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Abstract 


Tins  report  provides  a  technical  description  of  the  models  that  comprise  Component  1  of 
a  4-year  Chemical,  Biological,  Radiological  and  Nuclear  Research  and  Technology  Initia¬ 
tive  (CRTI)  Project  02-0093RD  entitled  “An  Advanced  Emergency  Response  System  for 
Chemical  Biological  Radiological  and  Nuclear  (CBRN)  Hazard  Prediction  and  Assessment 
for  the  Urban  Environment”  whose  primary  objective  is  the  development  of  an  advanced, 
fully  validated,  state-of-the-science  modeling  system  for  the  prediction  of  urban  flow  (i.e., 
turbulent  flow  through  cities)  and  the  concomitant  problem  of  modeling  the  dispersion  of 
CBRN  agents  released  in  a  populated  urban  environment.  The  principal  module  of  Com¬ 
ponent  1  is  urbanSTREAM,  which  is  a  general  second-order  accurate  finite-volume  code 
designed  for  the  simulation  of  urban  flow  using  a  two-equation  turbulence  closure  model 
(namely,  the  standard  k-e  model  and  the  limited-length-scale  k-(  model).  Component  1 
also  incorporates  a  module  (urbanGRID)  for  the  automatic  generation  of  grids  in  the  com¬ 
putational  domain  when  provided  with  detailed  geometric  information  on  the  shapes  and 
locations  of  buildings  in  the  urban  environment  in  the  form  of  Environmental  Systems 
Research  Institute  (ESRI)  Shapefiles.  Finally,  Component  1  also  includes  modules  for  the 
prediction  of  urban  dispersion  in  the  Eulerian  framework:  namely,  urbanEU  which  is  an  Eu- 
lerian  grid  dispersion  model  based  on  numerical  solution  of  a  A'-theory  advection-diffusion 
equation  (source-oriented  approach)  and  urbanAEU  which  is  a  receptor-oriented  dispersion 
model  based  on  numerical  solution  of  the  adjoint  of  the  A'-theory  advection-diffusion  equa¬ 
tion.  The  report  describes  the  current  status  of  the  urban  microscale  flow  and  dispersion 
modeling  system  developed  in  Component  1  and  presents  a  case  study  of  the  application 
of  this  system  in  stand-alone  mode  to  the  simulation  of  flow  and  tracer  dispersion  in  a 
real  cityscape.  The  numerical  simulations  of  Intensive  Observation  Period  (IOP)  9  in  the 
Joint  Urban  2003  (JU2003)  experiments  conducted  in  Oklahoma  City,  Oklahoma  provide 
an  initial  demonstration  that  the  developed  modeling  system  can  correctly  reproduce  many 
features  of  the  flow  and  dispersion  in  a  real  urban  environment. 
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1 


Resume 


Ce  rapport  fournit  une  description  technique  des  modeles  compris  dans  ia  Composante  l 
d’un  projet  de  4  ans  appele  Initiative  de  recherche  et  de  technologie  chimique,  hiologique.  ra- 
diologique  et  nuclaire  (IRTC),  Pro  jet  02-0093RD  intitule  “Systeme  aineliore  d'  intervention 
d’urgence  de  prediction  de  dangers  chimiques,  hiologiques.  radio logiques  et  nucleates  {CB- 
RN)  et  devaluation  du  milieu  urbain” .  L’objectif  principal  de  ce  projet  consiste  a  devel- 
opper  un  systeme  de  modelisation,  a  la  pointe  de  la  recherche  scientifique  et  totalement 
valide,  de  prediction  du  courant  urbain  {ex.  :  courant  turbulent  a  travers  les  villes)  ainsi 
que  du  probleme  concomitant  de  la  modelisation  de  la  dispersion  des  agents  CBRN  laches 
dans  le  milieu  urbain  peuple.  Le  module  principal  de  la  Composante  1  est  urbanSTREAM, 
un  code  genera]  de  deuxieme  ordre  de  volume  fini  exact  congu  pour  simuler  le  courant  ur¬ 
bain  en  utilisant  un  schema  de  fermeture  de  turbulence  a  deux  equations  (dont  le  schema 
jfc-e  et  le  schema  k-£  d’echelle  de  longueur  limitee).  La  composante  1  incorpore  aussi  un 
module  (urbanGRID)  pour  la  generation  automatique  de  grilles  dans  le  domaine  compu- 
tationnel  quand  sont  fournies  les  informations  geometriques  detaillees  des  formes  et  des 
emplacements  des  batiments  du  milieu  urbain  sous  forme  de  Fichiers  de  formes  provenant 
de  Tinstitut  Environmental  Systems  Research  Institute  (ESRI).  Enfin,  la  composante  1 
comprend  aussi  des  modules  de  prediction  de  la  dispersion  urbaine  dans  le  cadre  conceptuel 
eulerien,  dont  urbanEU  qui  est  un  modele  de  dispersion  de  grille  eulerienne,  base  sur  la 
solution  numerique  d'une  equation  advection-diffusion  de  la  theorie  K  (methode  axee  sur  la 
source)  et  urbanAEU  qui  est  un  modele  de  dispersion  axe  sur  Les  recepteurs,  base  sur  une 
solution  numerique  de  Tadjoint  de  l’equation  advection-diffusion  de  la  theorie  K.  Ce  rapport 
decrit  le  statut  actuel  du  systeme  de  courants  urbains  a  micro-echelle  et  de  modelisation 
de  la  dispersion  dans  la  Composante  1  et  presente  une  etude  de  cas  de  1'application  de  ce 
systeme,  en  mode  autonome,  a  la  simulation  des  courants  et  de  la  dispersion  des  traceurs 
dans  un  paysage  urbain  reel.  Les  simulations  numeriques  de  la  Periode  d 'observation  in¬ 
tensive  (POI)  9  des  experiences  Joint  Urban  2003  (JU2003),  conduites  a  Oklahoma  City 
en  Oklahoma,  commencent  a  demontrer  que  le  systeme  de  modelisation  mis  an  point  a 
la  capacite  de  reproduire  correctement  beaucoup  de  caracteristiques  des  courants  et  de 
dispersion  dans  un  milieu  urbain  reel. 
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Executive  summary 


Introduction:  The  release  of  chemical,  biological,  or  radiological  (CBR)  agents  by  ter¬ 
rorists  or  rogue  states  in  a  North  American  city  (densely  populated  urban  center)  and  the 
subsequent  exposure,  deposition,  and  contamination  are  emerging  threats  in  an  uncertain 
world.  As  a  consequence,  the  understanding  of  the  wind  flow  in  an  urban  area  and  the 
concomitant  dispersion  of  material  released  into  that  flow  is  crucially  important.  Unfortu¬ 
nately,  predictive  models  for  urban  flow  and  dispersion  do  not  exist.  In  view  of  this,  the 
work  described  in  this  report  was  undertaken  to  address  this  important  capability  gap. 

Results:  This  report  provides  a  technical  description  of  the  models  that  comprise  Compo¬ 
nent  1  of  a  4-vear  Chemical,  Biological,  Radiological  and  Nuclear  Research  and  Technology 
Initiative  (CRTI)  Project  02-G093RD.  This  project  entitled  “An  Advanced  Emergency  Re¬ 
sponse  System  for  Chemical  Biological  Radiological  and  Nuclear  (CBRN)  Hazard  Prediction 
and  Assessment  for  the  Urban  Environment”  involved  a  collaborative  model  development 
effort  by  Defence  R&D  Canada  -  Suffield  and  Enviroment  Canada.  The  primary  objective 
of  the  project  was  to  develop  an  advanced,  fully  validated,  state-of-the-science  modeling  sys¬ 
tem  for  the  prediction  of  urban  flow  and  for  the  modeling  of  the  dispersion  of  CBRN  agents 
released  in  these  flows.  Defence  R&D  Canada  -  Suffield  was  responsible  for  Component  1 
of  CRTI  Project  02-0093RD  which  involved  (1)  development  and  implementation  of  a  com¬ 
putational  fluid  dynamics  model  for  the  simulation  of  urban  flow  in  an  arbitrary  cityscape; 
(2)  development  of  a  grid  generation  capability  for  the  automatic  generation  of  grids  in  the 
computational  domain  when  provided  with  detailed  geometric  information  on  the  shapes 
and  locations  of  buildings;  (3)  provision  of  predictive  models  for  urban  dispersion;  and,  (4) 
validation  of  all  these  model  components.  The  report  provides  a  detailed  technical  descrip¬ 
tion  of  the  urban  flow  and  dispersion  models  developed  for  CRTI-02-0093RD  and  presents 
a  case  study  of  the  application  of  these  models  to  the  simulation  of  flow  and  contaminant 
dispersion  in  a  real  cityscape.  The  numerical  simulations  of  Intensive  Observation  Period 
(IOP)  9  in  the  Joint  Urban  2003  (JU2003)  experiments  conducted  in  Oklahoma  City,  Okla¬ 
homa  provide  an  initial  demonstration  that  the  developed  modeling  capability  can  correctly 
reproduce  many  features  of  the  flow  and  dispersion  in  a  real  urban  environment. 

Significance  and  Future  Plans:  The  availability  of  high-fidelity,  time-dependent  models 
for  the  prediction  of  a  CBR  agent’s  movement  and  fate  in  a  complex  urban  environment 
will  provide  the  strongest  technical  and  scientific  foundation  for  support  of  Canada's  more 
broadly  based  effort  at  advancing  counter-terrorism  planning  and  operational  capabilities. 
The  final  product  of  this  research  and  development  effort,  will  be  a  high-fidelity  multi¬ 
scale  CBR  modeling  system  that  will  be  fully  operational  at  the  Environmental  Emergency 
Response  Division  (EERD)  at  the  Canadian  Meteorological  Centre  (CMC).  This  resource 
will  serve  as  a  nation-wide  general  problem-solving  tool  for  first- res  ponders  involved  with 
CBR  incidents  in  the  urban  environment  and  beyond.  Future  plans  for  the  modeling  system 
described  here  include:  (1)  incorporation  of  more  advanced  turbulence  models  for  prediction 
of  complex  flows;  (2)  inclusion  of  thermal  effects;  and.  (3)  inclusion  of  an  adaptive  mesh 
refinement  capability  (among  others). 
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System:  Component  1  of  CRTI  Project  02-0093RD.  (DRDC  Suffield  TR  2007-067).  Defence 
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Sommaire 


Introduction:  La  liberation  d'agents  chimiques,  biologiques  ou  radiologiques  (CBR)  par 
des  terroristes  ou  des  etats  sans  scrupules,  dans  une  ville  nord-americaine  (centre  urbain  pe¬ 
ople  densement)  et  1 ‘exposition,  depot  et  contamination  atmospheriques  qui  s'en  snivraient. 
sont  les  menaces  emergentes  d’un  monde  incertain.  En  consequence,  la  comprehension 
des  courants  venteux  dans  une  zone  urbaine  et  de  la  dispersion  concomitante  des  matieres 
liberees  dans  ces  courants  est  tres  importante.  Les  modeles  de  prediction  des  courants  ur- 
bains  et  de  dispersion  n ’existent  malheureusement  pas.  De  ce  fait,  les  travaux  decrits  dans 
ce  rapport  ont  ete  entrepris  pour  aborder  le  probleme  de  cette  importante  lacune. 


Resultats:  Ce  rapport  fournit  une  description  technique  des  modeles  compris  dans  la  Com- 
posante  1  d’un  pro  jet  de  4  ans  appele  Initiative  de  recherche  et  de  technologic  chimique,  bi- 
ologique,  radiologique  et  nucleaire  (IRTC),  Pro  jet  02-0093RD.  Ce  pro  jet  intitule  “Systeme 
ameliore  d’intervention  d’urgence  de  prediction  de  dangers  chimiques,  biologiques,  radi¬ 
ologiques  et  nucleaires  (CBRN)  et  devaluation  du  milieu  urbain”  consiste  en  un  efTort  col- 
laboratif  de  developpement  de  modele  entre  R&D  pour  la  defense  Canada  -  Suffieid  et  Envi- 
ronnement  Canada.  L’objectif  principal  consistait  a  developper  un  systeme  de  model isat ion. 
a  la  pointe  de  la  recherche  scientifique  et  totalement  valide,  de  prediction  des  courants  ur- 
bains  et  a  modeliser  la  dispersion  d’agents  CBRN  liberes  dans  ces  courants.  R&D  pour  la 
defense  Canada  -  Suffieid  etait  charge  de  la  Composante  1  du  Projet  Project  02-0093RD 
de  IRCT  comprenant  (1)  le  developpement  et  ^implementation  du  modele  de  calcul  de 
la  dynamique  des  fluides  pour  la  simulation  des  courants  urbains  dans  un  paysage  urbain 
arbitraire;  (2)  le  developpement  d’une  capacite  a  generer  des  grilles  pour  la  generation  au- 
tomatique  de  grilles  dans  le  domaine  computationnel  quand  sont  fournis  les  renseignements 
geometriques  detail  les  des  formes  et  des  emplacements  des  batiments;  (3)  la  provision  de 
modeles  de  prediction  de  la  dispersion  urbaine;  (4)  la  validation  de  toutes  les  composantes 
de  ces  modeles.  Le  rapport  procure  une  description  technique  detaillee  des  courants  et 
des  modeles  de  dispersion  urbains  developpes  pour  IRTC-02-0093RD  et  presente  un  etude 
de  cas  pour  1’ application  de  ces  modeles  de  simulation  des  courants  et  de  la  dispersion 
des  contaminants  dans  un  paysage  urbain  reel.  Les  simulations  numeriques  de  la  Periode 
d’observation  intensive  (POI)  9  dans  les  experiences  Joint  Urban  2003  (JU2003)  conduites 
a  Oklahoma  City  en  Oklahoma,  commencent  a  demontrer  que  les  capacites  de  modelisation 
mises  au  point  sont  capables  de  reproduce  correctement  beaucoup  de  caracteristiques  des 
courants  et  de  dispersion  dans  un  milieu  urbain  reel. 


Portee  des  Resultats  et  les  Plans  Futurs:  La  dispombilite  de  modeles  de  haute  fidelite 
et  dependants  du  temps,  concernant  la  prediction  du  mouveinent  et  du  sort  des  agents  CBR 
dans  un  milieu  urbain  complexe,  procurera  les  fondetnents  techniques  et  scientifiques  les  plus 
sol  ides  en  soutien  aux  efforts  plus  generaux  du  Canada  a  ameliorer  la  plan  ificat  ion  du  con- 
treterrorisme  et  de  ses  capacites  operationnelles.  Le  produit  final  de  cet  effort  en  recherche 
et  developpement  consjstera  en  un  systeme  de  modelisation  CBR  a  multi  echelle  et  de  haute 
fideli’e  qui  sera  totalement  operationnel  a  la  Division  de  la  response  aux  urgences  environ- 
nementales  du  Centre  Meteorologique  Canadien  (CMC).  Cette  ressource  servira  d'outil  fie 
resolution  de  problemes  generaux.  a  1 ’echelle  de  la  nation,  aux  premiers  intervenants  par¬ 
ticipant  aux  incidents  CBR.  dans  le  milieu  urbain  et  au-dela.  Les  plans  futurs  du  systeme 
de  modelisation  d'ecrit  ici  comprend:  (1)  1’ incorporation  de  modeles  ameliores  de  turbu¬ 
lences  pour  la  prediction  des  courants  complexes;  (2)  l’inclusion  des  effets  thermiques;  et 
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(3)  1'inclusion  d’une  capacite  amelioreede  raffinement  des  mailles  (entre  autres). 
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Introduction 


Background 

The  environmental  and  toxiological  impact  of  the  downwind  transport  and  diffusion  of  con¬ 
taminants  released  into  the  atmosphere  has  become  increasingly  important  in  recent  years. 
Considerable  interest  has  been  focused  on  the  prediction  of  mean  concentration  levels  down¬ 
wind  of  contaminant  sources  in  the  turbulent  atmospheric  boundary  layer.  Consequently, 
atmospheric  transport  and  diffusion  models  have  played  an  important  role  in  emergency 
response  systems  to  toxic  releases  and  have  been  used  in  calculating  the  transport,  diffu¬ 
sion,  and  deposition  of  toxic  chemical,  biological,  or  radiological  materials  released  (either 
accidentally  or  deliberately)  into  the  turbulent  atmospheric  boundary  layer  over  relatively 
smooth  and  horizontally  homogeneous  surfaces.  For  example,  military  and  civilian  (gov¬ 
ernment  and  commercial)  emergency  response  models  commonly  use  standard  Gaussian 
plume  or  puff  models  which  employ  semi-empirical  relationships  for  plume  or  puff  growth 
with  the  mean  wind  and  turbulence  field  obtained  either  from  similarity  theory  or  from  the 
use  of  simple  diagnostic  wind  fields  constructed  from  interpolation  and/or  extrapolation 
of  sparse  observational  data.  The  advantages  of  these  approaches  for  wind  flow  specifica¬ 
tion  are  their  simplicity,  general  applicability  in  simple  atmospheric  conditions,  and  most 
importantly,  their  limited  computational  requirements.  While  this  approach  is  useful  for 
a  landscape  that  is  relatively  flat  and  unobstructed,  it  is  wholly  inadequate  for  surface- 
atmosphere  interactions  over  "complex”  surfaces  (viz.,  most  of  the  real  world)  such  as  cities 
and  other  built-up  areas. 

It  needs  to  be  emphasized  that  as  the  fraction  of  the  World’s  population  that  live  in  cities 
continues  to  grow,  it  is  becoming  increasingly  important  to  address  the  urgent  problem 
of  modeling  of  the  dispersion  of  toxic  releases  in  the  urban  environment,  characterized  by 
extremely  diverse  length  and  time  scales  and  complex  geometries  and  interfaces.  Indeed, 
a  typical  urban  canopy  consists  of  a  large  collection  of  buildings  and  other  obstacles  (e  g., 
cars  lining  a  street,  treed  areas  in  city  green  spaces,  etc.)  that  are  aggregated  into  complex 
structures.  When  this  rough  surface  interacts  with  the  atmospheric  flow  within  and  above 
it,  the  disturbed  flow  field  can  become  extremely  complex  (e.g.,  curved  mean  streamlines, 
large  velocity  gradients,  sharp  velocity  discontinuities,  flow  separations  and  reattachments, 
cavity  regions,  recirculation  zones,  and  strongly  inhomogeneous  turbulence).  Understanding 
the  complex  flow  of  the  wind  through  and  above  the  urban  environment  and  the  dispersion 
of  contaminants  released  into  that  flow  is  both  necessary  and  important.  In  view  of  this, 
we  require  physically- based  urban  wind  models  that  can  provide  the  needed  spatial  pattern 
of  urban  wind  statistics  required  to  "drive”  modeling  of  dispersion  of  contaminants  within 
the  street  canyons  of  an  urban  environment  (where  it  is  venting  of  these  street  canyons  that 
is  important  for  determination  of  the  contaminant  concentrations). 

Ttiis  identified  capability  gap  was  the  motivation  for  the  development  of  an  advanced  emer¬ 
gency  response  system  for  chemical,  biological,  radiological  and  nuclear  (CBRN)  hazard 
prediction  and  assessment  for  the  urban  environment  sponsored  by  Chemical,  Biological, 
Radiological  and  Nuclear  Research  and  Technology  Initiative  (CRTI)  under  Project  02- 
QG93RD.  The  principal  objective  of  this  project  was  to  develop  an  advanced,  fully  vali¬ 
dated,  state-of-the-science  modeling  system  for  the  predition  of  urban  flow  (he.,  turbulent 
flow  through  cities)  and  the  concomitant  problem  of  the  dispersion  of  CBRN  agents  re- 
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leased  into  these  complex  flows.  This  system  will  allow  the  dispersion  of  CBRN  materials 
to  be  modeled  over  a  vast  range  of  length  scales  at  the  appropriate  resolution  for  each 
scale:  namely,  in  the  near  field  (up  to  about  2  km)  where  dispersion  is  governed  by  the 
micro-scale  regime  of  the  planetary  boundary  layer:  to  the  intermediate  field  between  about 
2  and  20  km  where  dispersion  is  governed  by  the  local  or  tneso-7  scale:  through  the  far 
field  covering  the  range  from  about  20-200  km  (meso-/3  scale)  and  from  about  200-2000  km 
(meso-a  scale)  which  correspond  to  dispersion  at  the  regional  scale;  and.  finally  out  to  the 
very  far  field  encompassing  scales  greater  than  about  2000  km  corresponding  to  dispersion 
on  the  large  (synoptic  and  global)  scales.  The  development  of  this  proposed  high-fidelity, 
multi-scale  and  multi- physics  modeling  system  will  provide  a  real-time  modeling  and  simu¬ 
lation  tool  for  prediction  of  injuries,  casualities,  and  contamination  and  for  making  relevant 
decisions  (based  on  the  strongest  technical  and  scientific  foundations)  required  in  the  sup¬ 
port  of  Canada’s  more  broadly  based  efforts  at  advancing  counter-terrorism  planning  and 
operational  capabilities. 

The  multi-scale  modeling  system  for  emergency  response  consists  of  five  major  components 
shown  in  the  schematic  diagram  of  Figure  1.  These  five  components  can  be  described 
briefly  as  follows.  Component  l  involves  the  development  of  models  to  predict  the  mean 
flow  and  turbulence  in  the  urban  complex  at  the  microscale  (from  the  building  and  street 
scale  up  to  a  length  scale  of  about  2  km).  Two  kinds  of  models  have  been  developed  for  this 
purpose:  namely,  high-resolution  building-aware  models  for  urban  flow  where  buildings  are 
explicitly  resolved;  and,  virtual  building  models  for  urban  flow  where  groups  of  buildings 
are  represented  simply  in  terms  of  a  distributed  drag  force. 

Component  2  involves  the  inclusion  of  the  effects  of  urban  terrain  on  the  subgrid  scales 
of  a  mesoscale  meteorological  model  (Global  Environmental  Multi-scale  Local  Area  Model, 
or  GEM  LAM)  through  an  urban  parameterization.  This  parameterization  is  required  to 
account  properly  for  the  area-averaged  effects  of  form  drag,  increased  turbulence  production, 
heating  and  surface  energy  budget  modification  due  to  the  presence  of  buildings/ obstacles 
and  urban  land  use  within  the  urban  environment.  Component  3  involves  coupling  the  urban 
microscale  How  models  developed  in  Component  1  with  the  urbanized  mesoscale  model 
developed  in  Component  2.  The  interface  between  the  urban  microscale  flow  models  and 
the  urbanized  GEM  LAM  model  is  demanding  in  that  the  information  transfer  between  the 
two  models  must  honor  physical  conservation  laws,  mutually  satisfy  mathematical  boundary 
conditions,  and  preserve  numerical  accuracy,  even  though  the  corresponding  meshes  might 
differ  in  structure,  resolution,  and  discretization  methodology. 

Component  4  involves  using  the  mean  flow  and  turbulence  predicted  by  the  multi-scale 
flow  model  developed  in  Component  3  to  "drive’’  a  Lagrangian  Stochastic  (LS)  model  for 
the  prediction  of  urban  dispersion  of  CBRN  agents.  The  application  of  LS  models  to 
atmospheric  dispersion  in  general  (and,  urban  dispersion  in  particular)  is  recommended 
because  LS  models  (1)  are  (in  principle)  the  most  flexible  and  the  most  easily  able  to 
incorporate  all  the  known  statistical  details  on  the  complex  urban  flow  and  (2)  are  physically 
transparent,  and  easily  adapted  to  handle  particulates,  biological  or  radioactive  decay,  dry 
and  wet  deposition,  and  other  relevant  source  and  sink  mechanisms,  Finally.  Component  5 
involves  the  validation  of  the  multi-scale  modeling  system  for  both  the  urban  flow  and 
dispersion  components. 

A  series  of  reports  and  user’s  guides  describe  the  various  model  components  of  C'RTI  Project 
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02-0093 RD.  The  primary  objective  of  this  report  is  to  describe  the  present  status  (and,  more 
specifically,  the  technical  formulation)  of  the  ongoing  model  development  conducted  under 
Component  1  of  the  project.  In  addition,  the  predictive  capabilities  of  the  urban  microscale 
flow  and  dispersion  models  developed  in  Component  1  and  applied  in  stand-alone  mode  are 
demonstrated  by  presenting  results  from  a  case  study  involving  the  comparison  of  model 
predictions  with  experimental  data  obtained  from  a  comprehensive  full-scale  urban  field 
experiment  conducted  in  Oklahoma  City,  Oklahoma  in  July  2003  (Joint  Urban  2003,  or 
JU2003  field  experiment). 


Overview  of  the  urban  microscale  modeling  system 

The  urban  modeling  system  developed  for  Component  1  of  CRTI  Project  02-G093RD  in¬ 
cludes  five  main  modules:  urbanGRID,  urbanSTREAM,  urbanEU,  urbanAEU,  and  ur¬ 
ban  POST.  These  modules  and  how  they  interface  with  each  other  and  with  other  project 
components  are  shown  in  Figure  2, 

In  the  simplest  terms,  urbanGRID  imports  building  information  encoded  in  Environmental 
Systems  Research  Institute  (ESRI)  Shapefiles  and  uses  this  data  to  generate  a  structured 
grid  over  a  user  selected  computational  domain  in  a  given  cityscape.  Furthermore,  urban¬ 
GRID  imports  three-dimensional  meteorological  fields  (e.g.,  mean  wind,  turbulence  kinetic 
energy,  etc.)  provided  by  urban  GEM  LAM  and  uses  this  information  to  provide  inflow 
boundary  conditions  for  the  urban  microscale  flow  model. 

The  structured  grid  and  inflow  boundary  conditions  provided  by  urbanGRID  are  used  as 
input  by  urbanSTREAM  which  is  a  computational  fluid  dynamics  (CFD)  model  for  the 
numerical  simulation  of  the  flows  around  and  within  the  complex  geometries  of  buildings 
in  the  cityscape.  The  flow  solver  urbanSTREAM  provides  the  high-resolution  wind  and 
turbulence  fields  used  by  the  two  Eulerian  grid  dispersion  models  urbanEU  (source-oriented) 
and  urbanAEU  (receptor-oriented)  to  simulate  the  dispersion  of  contaminants  in  the  urban 
domain.  These  two  urban  dispersion  models  are  based  on  the  numerical  solution  of  K- 
theory  advectiomdiffusion  equation  or  its  adjoint.  Finally,  urbanPOST  is  used  to  process  the 
primary  output  files  from  urbanSTREAM  to  provide  an  appropriate  specification  of  wind 
statistics  required  as  input  by  either  the  two  Eulerian  urban  dispersion  models  urbanEU 
and  urbanAEU  or,  alternatively,  by  the  urban  Lagrangian  Stochastic  (LS)  particle  model 
urbanLS  (forward  dispersion  model)  or  urbanbLS  (backward  dispersion  model)  developed 
under  Component  4  of  CRTI  Project  G2-0093RD. 

The  memory  management  scheme  used  in  urbanSTREAM,  urbanEU,  and  urbanAEU  is  full 
dynamic  array  allocation  which  allows  the  maximum  array  dimensions  in  these  models  to 
be  easily  adjusted  to  match  the  requirements  of  a  particular  problem.  The  size  of  the  grid 
for  a  specific  application  is  determined  in  urbanGRID  which  provides  an  external  parameter 
file  that  specifies  the  required  array  sizes  for  all  the  major  arrays  in  urbanSTREAM.  The 
latter  model  then  dynamically  allocates  only  the  amount  of  memory  required  for  the  specific 
problem. 
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Technical  description  of  modules 


Grid  generator:  urbanGRID 

urbanGRID  is  a  program  that  is  concerned  with  the  automatic  generation  of  structured 
grids  (or  meshes)  in  the  computational  domain  for  numerical  solutions  involving  the  com¬ 
plex  geometry  of  buildings  and  other  obstacles  in  the  urban  environment  (real  cityscape). 
This  program  provides  a  predetermined  mesh  which  fills  the  entire  computational  domain, 
and  serves  as  the  input  data  for  spatial  discretization  required  to  “drive”  the  flow  solver 
program  urbanSTREAM.  The  mesh  determines  the  location  in  the  domain  at  which  the 
flow  quantities  will  be  evaluated  as  well  as  where  boundary  conditions  need  to  be  applied 
at  all  the  walls  and  roofs  of  arbitrary-shaped  buildings  (and  other  solid  surfaces)  in  the 
domain. 

The  geometric  information  on  shapes  and  locations  of  buildings  in  the  computational  domain 
(corresponding  to  a  region  of  a  real  city)  is  assumed  to  be  available  in  form  of  ESRI  Shape- 
files  [1].  An  ESRI  Shapefile  is  a  digital  vector  storage  format  that  stores  non-topological 
information  and  spatial  features  in  a  data  set,  with  the  geometry  For  each  feature  comprising 
a  set  of  vector  coordinates  and  associated  attribute  information*  A  Shapefile  is  actually  a 
set  of  files*  Three  individual  files  are  mandatory  and  these  store  the  core  data,  although 
in  general  a  Shapefile  may  also  include  up  to  eight  further  optional  individual  files*  urban¬ 
GRID  only  requires  Shapefile  information  stored  in  the  three  mandatory  files.  These  three 
individual  files  have  the  following  file  extensions  in  their  name:  (1)  *,shp  is  the  main  file 
and  is  a  direct-access,  variable  record-length  file  whose  records  store  the  feature  geometry 
(for  our  current  application,  the  key  features  are  area  features  corresponding  to  building 
“footprints"  encoded  as  closed-loop  double-digitized  polygons);  (2)  * .  shx  is  the  file  that 
stores  the  index  of  the  feature  geometry  and  contains  the  offset  of  the  main  file  record 
from  the  beginning  of  the  main  file;  and*  (3)  *  .dbf  is  a  file  in  the  form  of  a  dBASE  table 
that  contains  feature  attributes  (e.g.,  height  of  buildings)  with  one  record  per  feature*  there 
being  a  one-to-one  correspondence  between  the  geometry  and  attributes  based  on  record 
number. 

The  computational  domain  where  the  flow  field  will  be  calculated  is  specified  by  the  user 
who  provides  the  x  and  y  coordinates  (easting  and  northing  coordinates,  respectively)  of  the 
southwest  (x|^yert  ij^er)  and  northeast  y^£er)  corners  of  the  domain  in  the  Universal 

Transverse  Mercator  (UTM)  coordinate  system*  Additionally,  the  user  provides  the  easting 
and  northing  coordinates  of  the  southwest  (xg^er,  y^er)and  northeast  yfgj?er)  cor¬ 

ners  for  an  inner  domain  which  lies  strictly  inside  the  computational  domain  where  buildings 
are  explicitly  resolved  (building-aware  region)  in  the  sense  that  appropriate  boundary  con¬ 
ditions  are  imposed  on  all  the  building  surfaces  (e.g.,  wall,  roofs).  For  the  region  that  lies 
outside  the  inner  domain  but  inside  the  computational  domain,  the  buildings  are  treated 
as  virtual  in  the  sense  that  the  effects  of  these  unresolved  buildings  on  the  flow  are  rep¬ 
resented  simply  as  a  distributed  mean- momentum  sink  in  the  mean  momentum  equations 
(which  is  described  later  in  this  report).  urbanGRID  uses  the  Shapefile  C  Library  Version 
1.2  (Shapelib)  [2]  for  reading  ESRI  Shapefiles,  and  of  greatest  interest  here,  are  the  vector 
coordinates  of  all  the  vertices  (ordered  in  a  counter-clockwise  direction)  that  define  the 
footprints  of  the  buildings  (which  are  represented  in  UTM  coordinates)  and  the  attribute 
corresponding  to  the  heights  of  the  buildings  in  the  inner  domain*  Finally,  the  height  Zd  of 
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the  computational  domain  is  provided  as  a  user  input  to  urbanGRID. 


Once  the  computational  domain  has  been  defined,  a  structured  Cartesian  mesh  is  gener¬ 
ated  over  the  domain.  The  implicit  grid  structure  alleviates  the  need  to  store  the  mesh 
connectivity,  and  permits  the  use  of  rapid  iterative  solution  algorithms  for  the  computation 
of  the  flow,  which  makes  use  of  sparse  matrix  solvers,  urbanGRID  generates  a  Cartesian 
mesh  over  the  computational  domain  by  defining  an  arrangement  of  discrete  grid  points 
that  define  the  cell  faces  (grid  lines).  In  any  given  coordinate  direction,  the  spacing  of 
the  grid  points  starting  from  a  solid  surface  (e.g,,  ground  surface,  wall  or  roof  of  building, 
etc.)  is  stretched  as  one  moves  away  from  the  surface  along  the  coordinate  direction.  The 
stretching  function  used  in  urbanGRID  has  the  form 


(/3  +  2a)((0+l)/(^-l))l'-l,|-“,/(,^l-^  +  2a  .  „ 

"  (to  +  l)(l  +  ((/>  +  !)/(/>-  >/«*—>)  '  3  1,2 . '' 


(1) 


where  Lt  is  the  length  of  the  domain  in  the  i- th  coordinate  direction,  is  the  number  of 
grid  points  in  the  i- th  coordinate  direction,  x*[j[  is  the  j- th  grid  point  coordinate  in  the  i- th 
coordinate  direction  for  a  uniform  grid  spacing,  and  xSii[j]  is  the  corresponding  j-th  grid 
point  coordinate  in  the  i-th  direction  for  the  stretched  grid.  In  Equation  (1),  a  €  {0, 1} 
and  0  €  (Roc)  are  parameters  of  the  stretching  function  that  are  used  to  determine, 
respectively,  the  direction  (positive,  negative,  or  in  both  directions)  of  stretching  along  the 
coordinate  axis  and  the  magnitude  of  the  stretching.  More  specifically,  if  a  —  0  the  spacing 
between  the  grid  points  is  stretched  along  the  positive  direction  of  the  x*- axis,  if  a  =  1  the 
stretching  is  applied  along  the  negative  direction  of  the  x^-axis,  and  if  a  =  ^  the  stretching 
is  along  both  the  negative  and  positive  directions  of  the  xi-axis.  Finally,  the  stretching 
of  the  spacing  between  grid  points  along  a  coordinate  axis  increases  as  0  — *  1*,  whereas 
0  —  qq  results  in  a  uniform  spacing  between  grid  points. 


In  the  inner  region  of  the  computational  domain,  the  buildings  are  explicitly  resolved  and 
it  becomes  imperative  to  define  an  appropriate  data  structure  to  encode  information  about 
whether  a  control  volume  (grid  cell)  is  a  fluid  cell  (which  lies  mostly  or  completely  within 
the  fluid  region)  or  an  obstacle  cell  (which  lies  mostly  or  completely  inside  a  building), 
and  if  a  fluid  cell  whether  any  face  of  this  cell  abuts  against  a  wall  or  roof  of  a  building 
(the  latter  information  being  required  for  the  implementation  of  the  wall  function  boundary 
condition  in  urbanSTREAM).  To  this  purpose,  urbanGRID  uses  an  8- bit  integer  FLAG  array 
to  encode  this  information  as  follows: 


high  _ low 


C_D 

on/off 

CJJ 

on/ofF 

CJB 

on/off 

C_T 
on/ off 

c_s 

on/off 

CJI 

on/off 

C.W  C_E 

on/off  on /off 

1  ON 
0  OFF 

128 

64 

32 

16 

8 

4 

2  1 

If  bit  number  7  (C_0  bit)  of  the  FLAG  array  is  set  for  a  cell  with  indices  (viz., 

FLAG(i,j,fc)  =  C,G  =  64),  then  this  cell  is  an  obstacle  cell:  otherwise,  it  is  a  fluid  cell.  Bit 
number  8  (CJ>  bit)  of  the  FLAG  array  determines  whether  the  cell  lies  in  the  region  where 
buildings  are  considered  to  be  virtual  and,  hence,  represented  as  a  distributed  momentum 
drag  force  in  the  mean  momentum  equation.  More  specifically,  if  bit  8  is  switched  on  in 
the  FLAG  array  for  cell  then  this  cell  lies  within  the  “virtual  building'’  region  of 

the  computational  domain.  The  remaining  lower-order  bits  of  the  FLAG  array  determine 
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whether  a  fluid  cell  is  also  a  boundary  cell  in  the  sense  that  one  or  more  faces  of  this  cell 
borders  an  obstacle  cell  For  example,  a  Hu  id  cell  for  which  the  lowest  bit  (or.  C  J)  bit)  of 
the  FLAG  array  is  ON  signifies  that  the  east  face  of  the  corresponding  fluid  cell  abuts  against 
a  building  surface.  Similarly,  the  second  (C_W)?  third  (CJJ),  fourth  (CJS),  fifth  (C_T),  and 
sixth  (CJ3)  hits  of  the  FLAG  array  if  switched  ON  imply  that  the  west,  north,  south,  top,  or 
bottom  faces  of  the  associated  fluid  cell  borders  on  a  building  surface. 

Because  each  building  extracted  from  an  ESRI  Shapefile  has  a  constant  height  attribute, 
it  is  computationally  more  efficient  to  consider  the  intersection  of  the  projection  of  the 
cell  in  the  horizontal  x-y  plane  with  the  building  footprint  (represented  by  a  closed-loop 
polygon  in  the  Shapefile)  in  order  to  determine  whether  a  particular  grid  celt  is  a  fluid  cell 
or  an  obstacle  cell.  To  this  end.  a  bounding  box  (or,  smallest  rectangle  that  completely 
surrounds  the  building  footprint  and  has  sides  that  are  parallel  to  the  z-axis  and  ?/-axis)  is 
constructed  for  each  building  within  the  inner  (or,  building-aware)  region.  Clearly,  only  grid 
cells  that  intersect  any  of  these  bounding  boxes  in  the  inner  region  can  possibly  be  classified 
as  obstacle  cells  (with  the  appropriate  bit  in  the  associated  FLAG  array  set);  otherwise,  the 
grid  cell  must  necessarily  be  classified  as  a  fluid  cell. 

For  the  grid  cells  that  intersect  one  of  the  bounding  boxes  in  the  inner  region,  the  determi¬ 
nation  of  whether  these  cells  are  obstacle  or  fluid  cells  depends  on  whether  the  centroid  of 
the  cell  lies  inside  or  outside  the  building,  respectively.  A  ray-casting  approach  is  used  to 
ascertain  this  condition.  As  shown  in  Figure  3,  the  ray-casting  approach  for  determination 
of  whether  a  grid  cell  should  be  classified  as  a  fluid  cell  or  obstacle  cell  involves  casting  a 
ray  r  from  the  cell  centroid  c  and  simply  counting  the  number  of  intersections  of  r  with  the 
various  faces  of  the  building  enclosed  by  the  bounding  box  that  the  grid  cell  overlaps.  If  the 
grid  cell  centroid  lies  inside  a  building  (labeled  B  in  Figure  3),  the  number  of  intersections 
of  r  with  the  building  faces  must  be  odd.  Because  the  building  footprint  is  encoded  as  a 
list  of  vector  coordinates  of  vertices  forming  a  number  of  line  segments  that  constitute  a 
closed  polygon,  we  need  to  determine  if  there  is  an  intersection  of  ray  r  with  each  of  these 
line  segments  and  then  count  the  total  number  of  intersections. 

The  determination  of  whether  two  line  segments  intersect  is  based  on  the  computation  of  a 
signed  area.  To  this  purpose,  let  us  focus  on  whether  line  segment  ab  intersects  line  segment 
cd  as  depicted  in  Figure  4.  To  ascertain  this,  we  calculate  the  signed  areas  for  the  following 
four  triangles:  A(acd),  A  (bed),  A(acb)  and  A(adb).  The  signed  area  for  a  triangle  A(abc) 
whose  three  vertices  are  a,  b  and  c  (ordered  in  either  a  clockwise  or  counter-clockwise 
direction)  is  given  by 

(ax  ay  1  \ 

bx  by  1  |  ,  (2) 

cy  l  / 

where  a  —  (ar,av),  b  —  (bxi  by)  and  c  =  (crtcy).  With  reference  to  Figure  4,  the  two 
line  segments  ab  and  cd  intersect  each  other  if  and  only  if  the  following  two  conditions  are 
verified: 

Aaccj  x  <c  0, 

Aacb  ^  Aafjk  <r  0. 


While  the  test  for  intersection  of  two  line  segments  provided  by  Equation  (3)  is  straight¬ 
forward,  its  implementation  on  the  necessarily  finite  precision  of  the  hardware  on  a  digital 
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computer  is  not.  Problems  can  arise  when  a  zero-area  triangle  corresponding  to  the  case 
when  the  vertices  of  a  triangle  are  exactly  collinear.  Given  the  finite  precision  of  a  digi¬ 
tal  computer,  distinguishing  between  the  case  of  a  zero-area  triangle  occurring  because  of 
round-off  error  and  precision  problems  or  because  of  an  exact  col  linearity  is  difficult  and 
needs  to  be  handled  carefully.  To  deal  with  this  possibility,  urbanGRID  uses  an  adaptive 
exact  arithmetic  procedure  developed  by  Schewchuk  [3]  in  its  computation  of  signed  areas. 
The  procedure  used  in  urbanGRID  for  this  is  as  follows.  We  first  calculate  the  signed  area 
given  by  Equation  (2)  using  floating  point  arithmetic,  and  then  estimate  the  maximum 
possible  value  for  the  round-off  error  using  the  error  bound  derived  in  [3].  If  this  error  is 
larger  than  the  computed  signed  area,  then  the  area  is  re-calculated  using  adaptive  pre¬ 
cision  floating  point  arithmetic.  If  this  result  is  identically  zero  too,  then  we  resolve  the 
degeneracy  in  the  signed  area  computation  by  using  a  “tie- breaking”  algorithm  based  on  a 
virtual  perturbation  approach  advocated  by  Edelsbrunner  and  Muncke  [4],  The  basic  idea 
underpinning  this  approach  is  to  resolve  “ties”  in  the  vector  coordinates  of  the  triangular 
vertices  that  result  in  an  “exact”  zero  signed  area  by  breaking  this  degeneracy  through  the 
introduction  of  a  unique  and  ordered  perturbation  of  the  vertices.  This  results  in  a  non¬ 
degenerate  signed  area.  Since  the  perturbation  is  virtual,  no  geometric  data  is  modified. 
It  is  noted  that  this  refinement  is  merely  a  technicality  needed  to  resolve  possible  finite 
precision  collinearity  of  vertices  of  a  triangle.  The  virtual  perturbation  of  the  vertices  never 
affects  the  actual  coordinates  of  these  vertices  (which  are  always  held  to  finite  precision), 
but  nevertheless  enables  an  unambiguous  calculation  of  the  signed  area. 


The  distribution  of  flow  variables  needs  to  be  specified  (viz.,  Dirichlet  boundary  conditions 
need  to  be  defined)  at  the  in  flow'  boundary  planes  for  the  computational  domain  used  for 
flow  calculation  in  urbanSTREAM.  To  this  purpose,  urbanGRID  interpolates  the  gridded 
mean  velocity  and  turbulence  fields  provided  by  urban  GEM  LAM  (a  prognostic  mesoscale 
model  with  an  urban  parameterization  developed  by  Environment  Canada  for  CRTI  Project 
02-0093RD).  The  gridded  mean  velocity  and  turbulence  kinetic  energy  fields  provided  by 
urban  GEM  LAM  in  a  large  domain  that  includes  the  computational  domain  defined  in 
urbanGRID  are  linearly  interpolated  to  the  grid  nodes  in  the  inflow  boundary  planes  for 
urbanSTREAM.  This  corresponds  to  a  one-way  coupling  between  the  mesoscale  flow  model 
urban  GEM  LAM  and  the  microscale  (building-aware)  flow  model  urbanSTREAM. 


Flow  solver:  urbanSTREAM 


The  prediction  of  the  complex  flow  (mean  wind  and  turbulence  quantities)  through  and 
above  an  urban  canopy  consisting  of  groups  of  buildings  in  various  configurations  that  are 
representative  of  a  real  cityscape  is  based  on  the  Reynolds- averaged  Navier-Stokes  (RANS) 
equations  whereby  the  turbulent  flow  is  considered  as  consisting  of  two  components:  a  fluc¬ 
tuating  part  and  a  mean  or  average  part.  The  Reynolds-averaged  continuity  and  momentum 
equations  for  an  incompressible  adiabatic  (i.e.,  not  buoyancy-affected)  fluid  are  a  system 
of  partial  differential  equations  governing  mass  and  momentum  conservation  and  can  be 
expressed  in  Cartesian  coordinates  as  follows: 


Continuity 


du2 

dxt 


—  0, 


(4) 
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Momentum 


dut  diijiti 
dt  +  dxj 


dp  d2 * *Ui 
dxt  +  U  dx'j 
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~  ^-(u'u')  -  2eijknjUk, 


(5) 


where  the  Reynolds  averaging  of  a  quantity  is  denoted  by  drawing  a  bar  over  the  quantity. 
Summation  is  implied  by  repeated  indices.  Here,  ut  and  u\  are  the  mean  and  fluctuating 
velocities  in  the  Xt -direct ion,  respectively,  with  i  —  1,  2,  or  3  representing  the  west-east 
x,  south-north  y,  or  vertical  z  directions;  x*  =  (x.y,  z)  =  x;  t  is  time;  u,  =  (tt,  v,  w); 
u,  —  (u,  u,  u));  v  is  the  kinematic  viscosity  and  p  is  the  kinematic  mean  pressure  (with 
p '  used  to  denote  pressure  fluctuations).  The  last  term  on  the  right-hand-side  (RHS)  of 
Equation  (5)  is  the  Coriolis  acceleration  term  (caused  by  the  Coriolis  force)  where  is  the 
Earth’s  rotation  vector  and  is  the  permutation  symbol  (which  has  the  value  zero  if  any 
pair  of  subscripts  is  identical,  and  is  (-l)9  otherwise  where  q  is  the  number  of  subscript 
transpositions  required  to  bring  {ijk)  to  the  natural  order  (123)).  Rotation  is  defined  by 
the  right-hand  rule  where  positive  is  clockwise  when  looking  in  the  direction  of  Qj.  Note 
that  the  Cartesian  coordinate  system  fixed  on  the  Earth’s  surface  can  have  any  orientation 
with  respect  to 


Reynolds-averaging  the  Navier-Stokes  equation  gives  rise  to  the  so-called  kinematic  Rey¬ 
nolds  stresses  which  are  defined  as  the  tensor  u'pi'y  The  Reynolds  stresses  depend  on  the 
velocity  fluctuations  u'  ,  and  introduce  new  unknown  quantities  in  the  RANS  equations,  and 
therefore  these  equations  no  longer  constitute  a  closed  system.  In  order  to  close  the  system 
of  equations,  we  need  further  equations  describing  the  relationship  between  u[u'j  and  the 
quantities  u,  and  p  that  we  seek  to  determine.  This  is  known  as  the  closure  problem  in 
turbulence  modeling.  One  of  the  simplest  turbulence  models  for  u'  u'  involves  approximating 
the  Reynolds  stress  components  by  analogy  with  a  Newtonian  type  of  linear  constitutive 
relationship  between  the  turbulence  stress  and  the  mean  strain-rate  tensor.  This  model 
uses  the  Boussinesq  eddy- viscosity  approximation  (which  is  perhaps  the  simplest  coordinate 
invariant  relationship  between  stresses  and  strains)  and  is  given  by 


where  vt  is  the  kinematic  eddy  viscosity,  k  =  -u'u'  is  the  turbulence  kinetic  energy  (TKE), 
and  Sl}  is  the  Kronecker  delta  function.  It  is  interesting  to  note  that  with  the  adoption  of 
the  Boussinesq  eddy-viscosity  approximation,  the  RANS  equations  have  the  same  form  as 
the  Navier-Stokes  equations  except  that  the  kinematic  molecular  viscosity  u  is  replaced  by 
an  effective  kinematic  viscosity  i/eff  =  (v  +  vt)  and  provided  the  mean  pressure  is  interpreted 
to  be  a  modified  mean  pressure1 . 


With  the  Boussinesq  eddy-viscosity  model  (EVM)  for  the  Reynolds  stresses,  the  task  of 
turbulence  closure  reduces  to  the  determination  of  the  eddy  viscosity.  Dimensional  analysis 
dictates  that  the  eddy  viscosity  i>t  be  determined  by  the  product  of  a  turbulence  velocity 
scale  uq  and  a  turbulence  length  scale  /q,  so 


=  iVo, 

1  When  Equation  (6)  is  inserted  into  Equation  (5),  the  term  \kStJ  in  the  Boussinesq  ap¬ 

proximation  for  the  Reynolds  stresses  is  usually  incorporated  with  the  pressure  p  to  define 

a  modified  pressure  p  +  (so,  p  — *  p  +  |/i:  =  pmod  in  the  RANS  equations). 
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where  u0  and  l0  can  vary  significantly  in  space  and  time  for  a  given  turbulent  flow  (and,  in 
particular,  for  an  urban  flow).  In  the  most  common  approach,  the  turbulence  velocity  and 
length  scales  are  based  on  the  turbulence  kinetic  energy  k  and  its  dissipation  rate  e  (i.e., 
rate  at  which  TKE  is  converted  into  thermal  internal  energy),  with  the  turbulence  velocity 
and  length  scales  chosen  to  be 


uo  =  cyv",  I,  =  (7a) 

Hence,  in  this  so-called  k-e  modeling  framework  at  high-Reynolds  (high-Re)  numbers,  the 
eddy  viscosity  is  given  by  (on  combining  turbulence  velocity  and  length  scales) 

k2 

=  >  (76) 

where  C M  is  a  closure  constant. 


The  modeled  transport  equation  for  the  turbulence  kinetic  energy  k  is 


where  Pk  is  the  production  of  k  defined  as 


Pk  = 


—dux 

ljdx  / 


(9) 


The  terms  in  the  modeled  TKE  transport  equation  [Equation  (8)]  are  from  left  to  right: 
local  rate  of  change  of  k,  transport  of  k  by  mean  advection,  transport  of  k  by  viscous 
and  turbulent  diffusion,  rate  of  production  of  k ,  and  rate  of  destruction  of  k.  Within  the 
framework  of  a  two-equation  turbulence  closure  scheme,  the  turbulent  diffusion  term  in  the 
exact  transport  equation  for  TKE  is  modeled  using  a  gradient  diffusion  hypothesis  as 


-u'u'u'  +  p'u\6ij 


i/t  dk 

ak  dij  ’ 


(10) 


where  ak  is  a  constant  that  functions  like  an  effective  Schmidt  number  for  the  turbulent 
diffusion  of  k.  Note  that,  with  this  model,  both  the  turbulent  and  molecular  diffusion  terms 
can  be  combined  into  a  gradient  transport  model  with  effective  viscosity  (is  +  ist/ak).  This 
resulting  simple  form  of  the  modeled  transport  equation  for  k  is  the  obvious  appeal  of  this 
formulation. 


An  exact  transport  equation  for  e  can  be  written,  but  this  equation  is  not  useful  owing  to 
the  fact  that  most  of  the  terms  in  the  equation  are  not  in  closed  form.  In  its  most  usual 
form,  the  modeled  equation  for  the  TKE  dissipation  rate  e  is 


dt  duj€ 
dt  +  dij 


+  -^{C(\Pk  -  C(2e), 


(11) 


where  Ct\,  Ct2  and  a(  are  closure  constants  to  be  determined, 
turbulent  diffusion  of  dissipation  (which  has  been  grouped  with 
term)  is  given  by 

Dl  =  —  \—  — 

*  dij  at  dij 


In  Equation  (11),  the 
the  molecular  diffusion 


(12) 
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where  af  is  a  constant  that  acts  like  an  effective  Schmidt  number  for  the  turbulent  diffu¬ 
sion  of  TKE  dissipation.  In  Equation  (11),  the  term  on  the  RHS  that  involves  C€\  is  the 
production  of  dissipation  term,  whereas  the  term  involving  C€ 2  is  the  destruction  of  dissi¬ 
pation  term.  The  factor  efk  in  the  production  and  destruction  terms  makes  these  terms 
dimensionally  correct  in  the  transport  equation  for  e. 


Within  the  Boussinesq  type  of  eddy- viscosity  approximation,  the  Reynolds  stresses  appear¬ 
ing  in  the  production  Pk  [cf.  Equation  (9)]  are  modeled  in  accordance  with  Equation  (6) 
resulting  in  the  following  modeled  form  for  P^: 


where 


/  dUi  dUj  \  dUi  _ 
k  Vl  \5zj  dxj )  dxj  " 

|5|  = -(25ii5ii)1/2,  = 


=  C„e\S\2, 

(13) 

(  9uj  \ 

\dxj  dxj' 

(14) 

Here,  StJ  is  the  mean  rate-of-strain  tensor  (or,  equivalently,  the  symmetric  part  of  the  mean 
velocity  gradient  tensor). 


Together,  the  transport  equations  for  fc  and  e  contain  five  closure  constants  (CMt  a cr€, 
C£  and  Ca)  which  must  be  determined  before  the  equations  can  be  solved.  Two  different 
high-Re  number  k-e  models  are  considered  here.  The  first  is  the  standard  k-e  model  and 
following  the  recommendation  of  Launder  and  Spalding  [5],  the  closure  coefficients  for  this 
turbulence  model  are 

C v  =  0.09,  ak  =  1,  -  T3,  Ccl  =  1.44,  Ce2  =  1.92.  (15) 

Here,  2  is  determined  from  the  experiment  on  the  decay  of  isotropic  turbulence,  and  C€y 
is  determined  experimentally  based  on  the  local  equilibrium  wall  shear  flow.  The  value  of 
ut  is  chosen  to  ensure  the  correct  log-law  behaviour  in  boundary-layer  flows.  The  second  k-e 
model  implemented  is  the  limited- length-scale  k-e  model  proposed  by  Apsley  and  Castro 
[6].  In  this  modified  k-e  model  the  length-scale  determining  equation  (i.e.,  e-transport 
equation)  is  modified  in  order  to  acknowledge  the  upper  bounds  on  the  turbulence  length 
scale  imposed  by  the  finite  depth  of  the  atmospheric  boundary  layer.  The  modification  is 
both  simple  and  elegant;  namely,  the  coefficient  Cel  in  the  production  of  dissipation  term 
is  replaced  by 

C£l^C£l  +  {Ca-C£l)-^-,  (16) 

^O.msix 

where  /o,max  —  h/3  is  the  maximum  turbulence  length  scale  imposed  by  the  depth  h  of  the 
atmospheric  boundary  layer  [and  Iq  is  defined  in  Equation  (7a)].  It  is  noted  that  with  the 
exception  of  Jo, max  1  the  modification  introduces  no  new  closure  constants  into  the  model 
(viz.,  the  closure  constants  C^r  tr*,  C€\.  and  C£ 2  assume  exactly  the  same  values  as  in 
the  standard  k-e  model). 


A  closed-form  solution  for  the  k-e  model  can  be  obtained  for  the  neutral  wall  shear  layer. 
The  solution  gives 


u* 


u  —  log  £  4-  k  — 


uz 


£  = 


kvz' 


(17) 
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where  u.  =  (— u'u/ )  is  the  friction  velocity  and  B  is  a  constant  of  integration.  For  this 
analytical  solution,  we  find  an  implied  value  for  the  von  Karman  constant,  kv,  of 

kv  —  x/Cj,  (C<2  C«i  (18) 

Using  the  closure  coefficient  values  for  the  fc-e  model,  kv  assumes  a  value  of  0.43.  The 
experimental  values  for  kv  are  primarily  in  the  range  0.41  ±  0.2,  so  the  implied  value  of  kv 
in  the  k-e  model  is  consistent  with  these  measurements.  The  stress-intensity  ratio,  u2.Jk,  is 
predicted  to  be  0.3  (using  Ch  =  0.09)  which  agrees  well  with  experimental  measurements 
of  this  quantity  in  many  shear  flows  [7]. 

The  technical  details  of  urbanSTREAM  used  in  the  high- resolution  computational  fluid 
dynamics  simulation  of  urban  flow  where  all  the  buildings  in  the  computational  domain 
are  explicitly  resolved  in  the  sense  that  appropriate  boundary  conditions  are  imposed  at  all 
the  building  surfaces  (e.g.,  walls,  roofs),  have  been  described.  However,  the  computational 
demands  of  such  an  approach  are  high  and,  as  a  consequence,  sufficiently  prohibitive  as 
to  preclude  the  simulation  of  flow  through  and  above  all  buildings  over  a  relatively  large 
computational  domain  in  a  typical  city  environment. 

In  lieu  of  imposing  correct  boundary  conditions  on  the  true  (but  usually  complex,  three- 
dimensional)  geometry  of  the  building  surfaces  and  fully  resolving  the  detailed  intricate  flow 
around  every  individual  building  in  a  computational  domain,  it  is  convenient  to  consider 
the  prediction  of  statistics  of  the  mean  wind  and  turbulence  in  an  urban  canopy  that  are 
obtained  by  averaging  horizontally  the  flow  properties  over  an  area  that  is  larger  than  the 
spacings  between  the  individual  roughness  elements  comprising  the  urban  canopy,  but  less 
than  the  length  scale  over  which  the  roughness  element  density  changes.  Towards  this 
objective,  we  consider  the  use  of  a  virtual  building  concept  in  urbanSTREAM  whereby 
groups  of  buildings  are  unresolved  and  the  aggregate  of  these  individual  buildings  is  treated 
simply  as  a  porous  medium.  Here,  the  effects  of  the  unresolved  buildings  {or,  virtual 
buildings)  on  the  flow  are  represented  through  a  distributed  mean- momentum  sink  in  the 
mean  momentum  equation.  Furthermore,  the  additional  source/sink  terms  in  the  supporting 
transport  equations  for  k  and  e  are  included  in  order  to  model  the  effects  of  the  virtual 
buildings  on  the  turbulence. 


The  foundations  of  a  systematic  mathematical  formulation  for  the  derivation  of  the  gov¬ 
erning  equations  for  flow  through  an  urban  canopy  where  the  individual  buildings  are  un¬ 
resolved  (virtual  buildings)  and  the  aggregate  of  various  groups  of  buildings  (and  other 
obstacles)  is  treated  simply  as  a  porous  medium  is  described  in  Lien,  Yee  and  Wilson  [8]. 
The  representation  of  groups  of  buildings  as  a  porous  medium  requires  an  additional  mod¬ 
eling  effort.  The  transport  equation  for  the  mean  momentum  is  modified  with  the  inclusion 
of  a  mean  drag  force  fd.i  due  to  the  effects  of  the  virtual  buildings  on  the  flow,  so 


dui  duiUt  dp  32u,  d  .-y— .-.  „  _  - 

aT  +  -ai7  =  -ai;+,'aJf  -  ~  +  fd-‘’ 

where  /d.,  is  parameterized  as  follows: 

fd.i  ~  -CdA  -f  , 


(19) 


(20) 
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where  Cp  is  the  drag  coefficient,  A  is  the  frontal  area  density  (frontal  area  of  buildings 

I/O 

exposed  to  the  wind  per  unit  volume),  and  Q  =  (u,Uj)  “  is  the  magnitude  of  the  (spatially- 
averaged)  time- mean  wind  speed.  Utilizing  the  eddy  viscosity  model  for  turbulent  stresses 
of  Equation  (6)  in  Equation  (20)  yields  the  following  final  form  for  the  mean  drag  force 
imposed  by  the  virtual  buildings  on  the  flow: 


fd.i 


(21) 


Standard  turbulence  models  need  additional  modifications  when  applied  to  urban  flows 
where  buildings  are  treated  as  virtual  (and  represented  as  a  distributed  mean-momentum 
sink).  Additional  modifications  required  in  the  transport  equations  for  k  and  e  were  de¬ 
rived  in  [8)  and  are  simply  summarized  here.  The  k-  and  e-equations  require  additional 
source/sink  terms  and  now  take  the  following  form: 


dk  dujk 
dt  +  dxj 


+  (Pfc  +  F)  -  e, 


(22) 


and 


de  duj£ 
dt  dr. 


dxj 


J?+£)^M(c*'(p‘+F>-c-4  <23> 


The  term  F  in  Equation  (22)  can  be  interpreted  as  an  additional  physical  mechanism  for 
the  production/dissipation  of  k ,  and  as  shown  in  [8]  can  be  approximated  as 


F  =  -CdA 


2 Qk  +  ~  ( tiiUfcu'ui.  )  +  ~  (  Ufett'uX  ) 


(24) 


where  the  triple  correlation  term  u'u'u't  here  is  modelled,  following  Daly  and  Harlow  [9],  as 


UiUtUh  ~  2C3  — 
e 


dk 


du'-ui 


- 1-  u'u't  -  *  ■ 

k  ldxi  1  1  3xt 


(25) 


where  the  closure  constant  Cs  ^0.3  is  used  in  the  present  study.  The  double  correlation 
(turbulent  stress)  u\uf}  in  Equations  (24)  and  (25)  are  modelled  using  the  Boussinesq  stress- 
strain  relationship  in  Equation  (6). 


Urban  dispersion:  urbanEU  and  urbanAEU 

Prediction  of  the  wind  field  statistics  using  urbanSTREAM  is  in  principle  pre-requisite  to 
or  co-requisite  with  prediction  of  the  simpler  (but  nevertheless  complex)  problem  of  scalar 
(contaminant)  diffusion  within  the  urban  environment.  In  particular,  the  prediction  of 
turbulent  flow  statistics  through  a  built-up  area  is  required  to  provide  the  input  information 
needed  to  "drive1'  a  physically- based  model  for  urban  dispersion.  Within  the  context  of 
Component  1  of  CRTI  Project  02-0Q93RD,  two  different  (but  complementary)  models  for 
urban  dispersion  from  the  Eulerian  perspective  have  been  developed,  and  will  be  described 
in  this  sub-section. 
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The  first  urban  dispersion  model,  referred  to  as  urbanEU,  uses  a  source-oriented  dispersion 
modeling  technique.  The  model  estimates  dispersion  about  the  Row  within  the  urban  canopy 
(viz.,  at  the  street  and  neighborhood  scale)  for  passive  scalars  released  into  that  flow.  In 
consequence,  the  mean  momentum  equation  is  decoupled  from  the  scalar  transport  equation 
with  the  velocity  field  independent  of  the  scalar  concentration  field.  More  specifically,  the 
additional  transport  equation  required  for  the  conservation  of  the  scalar  {the  contaminant) 
has  the  following  form: 


dC  dtijC 
dt  +  dxj 


{K  +  Kt) 


dc 

dxj 


+  Q* 


(26) 


where  C  is  the  concentration  of  a  contaminant  (g  m“3),  K  is  the  molecular  kinematic 
diffusivity  of  the  contaminant  in  air  (m2  s_1),  Kt  is  the  turbulent  (or,  eddy)  kinematic 
diffusivity  (m2  s^1),  and  Q  is  the  source  density  distribution  for  the  contaminant  (g  s-1 
m_3)« 


In  Equation  (26),  a  gradient  diffusion  model  is  used  as  perhaps  the  simplest  closure  model 
for  the  turbulent  scalar  fluxes,  although  this  assumption  has  not  been  rigorously  justified  for 
urban  dispersion.  In  particular,  within  this  popular  framework  the  turbulent  scalar  fluxes 
are  modeled  as 


u'c'  -  -Kt 


dC 

dxj  7 


(27) 


and  the  turbulent  diffusivity  Kt  is  obtained  from  the  turbulent  viscosity  i/t  (predicted  by 
urbanSTREAM)  in  combination  with  a  turbulent  Schmidt  number  Sct  in  the  following 
manner; 


(28) 


In  urbanEU,  we  adopt  this  concept  using  a  constant  turbulent  Schmidt  number  Sct  with 
a  value  of  0.63,  This  value  for  Sct  was  chosen  with  reference  to  Project  Prairie  Grass, 
a  benchmark  field  experiment  documenting  turbulent  dispersion  from  a  continuous  point 
source  near  the  ground  into  the  atmospheric  surface  layer,  where  it  has  subsequently  been 
shown  that  an  Eulerian  dispersion  model  is  in  excellent  conformance  with  these  high-quality 
experimental  measurements  provided  the  turbulent  Schmidt  number  has  a  value  3c*  ^  0.63 
{Wilson  [10]). 


In  the  source-oriented  dispersion  model  implemented  by  urbanEU,  the  goal  of  the  dispersion 
modeling  is  to  calculate  the  mean  concentration  "seen*  by  a  detector  at  a  given  receptor 
location  when  provided  the  source  distribution  (source  term)  Q  of  the  contaminant.  In 
other  words,  in  the  source-oriented  approach  the  advection-diffusion  equation  given  by 
Equation  (26)  is  solved  forward  in  time  for  a  given  source  distribution  of  contaminant  Q(x ,  t) 
to  obtain  the  mean  concentration  field  C(x,t).  Then,  the  concentration  of  contaminant 
'seen'  by  a  detector  at  a  given  receptor,  denoted  <J>(C),  can  be  defined  as  the  following 
integral  of  concentration  C(x,  t)  filtered  in  space  and  time  by  the  spatial-temporal  filtering 
function  h(x,t)  corresponding  to  the  response  function  of  the  detector: 

$(C)  =  j  C(x,  f)/i(x,  t)  dx  dt  =  { C.h ).  (29) 

Therefore,  this  integral  is  simply  a  linear  functional  of  C  and  expresses  averaging  of  the 
concentration  field  over  the  spatial  volume  and  averaging  time  imposed  by  the  detector  at. 
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the  receptor  location.  This  linear  functional  of  C  is  defined  by  the  inner  product  (C,  h) 
given  in  Equation  (29),  and  allows  one  to  calculate  detector  concentrations  at  any  number 
of  receptors  located  within  the  computational  domain.  However,  for  any  new  release  sce¬ 
nario  involving  a  different  source  distribution  Q.  the  solution  of  Equation  (26)  needs  to  be 
calculated  again  in  order  to  determine  $. 


Alternatively,  when  the  detector  concentration  at  a  fixed  receptor  is  of  principal  interest  for 
a  range  of  different,  emission  scenarios,  the  alternative  receptor-oriented  dispersion  modeling 
approach  is  useful.  As  a  result,  we  have  implemented  a  receptor-oriented  dispersion  model  in 
urbanAEU,  which  solves  the  following  adjoint  advect  ion-diffusion  equation  for  the  influence 
(or,  adjunct)  function  C '  (m-3): 


dC'  duj  C’ 
dt  dxj 


d 

dxj 


{K  +  I<t) 


+  h  , 


(30) 


where  h  is  the  space- time  filtering  function  of  the  detector  at  the  fixed  receptor  location. 
It  is  important  to  note  that  Equation  (30)  is  the  adjoint  of  Equation  (26)  with  C  and  C* 
related  to  each  other  through  the  source  function  Q  and  detector  (receptor)  function  h  as 
follows: 

<f>(C)  =  J  C*{x,t)Q(x,t)dxdt=  (C\Q)  =  (C,h).  (31) 

Equation  (31)  is  a  duality  relation  that  relates  C  with  C*  through  Q  and  h  and,  in  this 
sense,  C  and  C*  can  be  interpreted  as  dual  quantities.  In  the  receptor-oriented  approach, 
Equation  (30)  is  solved  backwards  in  time  for  a  given  receptor,  and  the  concentration  'seen' 
by  a  detector  at  this  receptor  can  now  be  calculated  directly  using  Equation  (31)  once  the 
source  distribution  Q  has  been  specified.  It  should  be  noted  that  these  computations  can 
be  repeated  for  any  source  distribution  Q  to  obtain  the  concentration  at  the  given  receptor 
without  having  to  re-solve  Equation  (30)  for  C*.  However,  a  new  adjunct  function  C*  needs 
to  be  determined  if  the  concentration  is  desired  for  a  different  receptor. 


Numerical  implementation 


Finite-volume  discretization 


The  coupled,  nonlinear  system  of  six  partial  differential  equations  which  models  unsteady, 
three-dimensional  turbulent  flow  field  in  urbanSTREAM  and  the  partial  differential  equa¬ 
tion  describing  scalar  transport,  in  urbanEU  and  urbanAEU  are  solved  using  a  finite- volume 
method.  All  these  partial  differential  equations  have  the  following  generic  form: 


d0 

dt 


du}0 

dxj 


(32) 


for  a  scalar  <p  (which  can  be  identified  with  any  of  the  mean  Cartesian  velocity  components 
ilj.  turbulence  kinetic  energy  k.  viscous  dissipation  e,  mean  concentration  C,  or  influence 
function  C*);  T*  is  the  diffusi vity  and  S#  represents  the  source/sink  term.  If  the  current 
time  level  for  0  is  denoted  by  the  superscript  n  +  1  and  the  time  increment  by  At.  a  two- 
or  three-level  time  discretization  can  be  used  to  approximate  the  local  rate  of  change  of  0 
as  follows: 

(do  \  1 

'dt)  ^  2At  {a°l^rl+I  ai<t>T'  +  ai<t>n~X )  ■ 
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where  a0  —  -a!  =  2,  a2  =  0  for  first-order  accuracy  (using  only  two  time  levels)  and  a0  -  3, 
a,  —  —4,  02  =  1  for  second-order  accuracy.  Both  options  are  available  in  urbanSTREAM 
and  urban(A)EU  for  approximating  temporal  derivatives.  Finally,  in  the  discretization  to 
follow,  it  is  assumed  that  the  source/sink  term  3$,  which  may  be  a  nonlinear  function  of  0, 
is  linearized  as  follows: 

5^  ~  Sp4>  +  Sc,  (34) 

where  it  is  implicitly  assumed  that  Sp  is  chosen  to  be  unconditionally  negative. 

Integration  of  the  transport  equation  governing  the  flow  property  0  in  Equation  (32)  over 
the  finite  volume  (or,  control  volume)  shown  in  Figure  5,  use  of  the  time  discretization  of 
Equation  (33)  to  approximate  the  local  tendency  term,  and  application  of  the  divergence 
theorem  results  in  a  discretized  transport  equation  for  0  that  expresses  a  balance  of  local 
flow  property  change,  advective  and  diffusive  fluxes  of  the  flow  property  across  the  cell  faces, 
and  generation/ destruction  of  0  embodied  by  the  volume- integrated  net  source/sink  term. 
The  resulting  full  implicit  time  discrete  version  of  Equation  (32)  takes  the  following  form 
[on  insertion  of  the  linearization  of  S$  given  by  Equation  (34)]: 

aPrP+l  =  Y  “nbflC1  +  Qc,  (35a) 

nb 

where 

Qc  =  ScAfl  --  +  020"  ']  (356) 

(Afl  is  the  volume  of  the  cell  with  center  node  P);  nb  refers  to  the  six  neighbors  of  cel! 
center  P  associated  with  nodes  at  the  cell  centers  to  the  west  (W),  east  (E),  south  (5), 
north  ( N ),  bottom  ( B ),  and  top  ( T )  of  the  center  node  P;  and,  the  six  surfaces  (faces)  of  the 
control  volume  centered  at  P  will  be  denoted  using  lower-case  letters  corresponding  to  their 
direction  { w ,  e,  s,  n,  b  and  t)  with  respect  to  P  (see  Figure  5).  The  surface  area  of  a  cell 
face  in  direction  k  relative  to  the  cell  center  P  will  be  denoted  by  A*  (k6  {u/s  e,  s,  n,  6,  t}). 


The  influence  coefficients  (not  to  be  confused  with  the  influence  function  C*)  an b  (nb  e 
{W\  E ,  S,  N,  B.T})  of  the  neighboring  unknowns  of  the  flow  property  0  take  the  following 
form,  assuming  that  the  diffusive  fluxes  have  been  approximated  by  central  differences  and 
the  advective  fluxes  have  been  approximated  using  a  first-order  upwind  scheme: 


aw 

a,£ 

“S 

aw 

aa 


I"1  fh  .  n i  *i 


yt/  p 

SXp£ 


w  +  max(uwAw,0), 
+  max(— ti(.Ae,0), 

F<ft.jjAs 

sVsp  . 

F  fjjtn  A„ 

.  ^Vpw  - 

r  Q'b 

&zWp 
A; 


+  max(53A„0), 

+  max(— t/nAn,0), 
+  max(u>t,A{,,  0), 

+  inax(-i5(At,  0), 


(36) 
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where  Sx\yp  =  xp  —  xw  denotes  the  distance  along  the  x-coordinate  direction  from  the 
west  neighbor  node  W  to  the  cell  center  node  P,  with  analogous  definitions  for  other  similar 
quantities  in  Equation  (36).  Finally,  the  central  coefficient  aP  is 

EAnAil 

anb  +  *2^ - SpAQ.  (37) 


urbanSTREAM  employs  a  fully-collocated,  cell-centered  storage  arrangement  for  all  trans¬ 
ported  flow  properties  as  well  as  for  pressure.  Furthermore,  urbanSTREAM  allows  two 
user-selected  options  for  the  approximation  of  the  advective  fluxes  across  a  cell  face;  namely, 
the  upwind  scheme  summarized  above  and  the  Upstream  Monotonic  Interpolation  for  Scalar 
Transport  (UMIST)  scheme  [11).  The  UMIST  scheme  is  a  total  variation  diminishing 
(TVD)  variant  of  the  higher-order  quadratic  upwind  interpolation  for  convective  kinematics 
(QUICK)  scheme  developed  originally  by  Leonard  [12]  to  approximate  volume-face  fluxes. 
This  second-order  scheme  for  the  discretization  of  the  net  advective  flux  through  a  con¬ 
trol  volume  combines  the  second-order  accuracy  of  a  central  differencing  scheme  with  the 
stability  inherent  in  an  upwind  differencing  scheme  by  using  in  each  direction  separately  a 
quadratic  upstream  interpolation.  The  second-order  accuracy  of  the  QUICK  scheme  mini¬ 
mizes  numerical  diffusion  errors  that  are  characteristic  of  first-order  accurate  discretization 
schemes  such  as  the  upwind  scheme.  The  UMIST  scheme  generalizes  QUICK  by  intro¬ 
ducing  a  limiter  to  give  a  monotonic  version  of  QUICK  that  satisfies  TVD  contraints.  In 
urbanSTREAM,  the  higher-order  advection  scheme  UMIST  is  implemented  simply  as  a 
first-order  upwind  approximation  (described  above)  augmented  with  deferred  corrections 
that  constitute  additional  contributions  to  the  source  term  Qc  given  in  Equation  (35b). 


To  this  purpose,  an  additional  source  contribution  (deferred  corrector  ‘source’)  Q^c  is 
added  to  Qc  (Qc  —*  Qc  +  QBC)  *n  ^le  implementation  of  UMIST.  This  deferred  corrector 
source  has  the  following  explicit  form: 


Qc°  =  0.5 1  [(u)*tf{r+)  -  (u)e  d(rc  )]  (<j>E  -  0P) 

~  [(“^(O  -  007^)]  (<f>p  -  M 

+  [(^)Tt^(rn)  ~  (")~l3(rn  )]  (^JV  -  4>p) 

-  [(*Os+tf«)  -  (®)7*(oj  typ  -  <t>s) 

+  [(«>)(+i?(r+)  -  (tB)'tf(rr)]  {<j>r  ~  Qp) 

-  [Hfctf(r^)  -  (w)“i?(r^)](^P  (38a) 

where  t?(r)  is  defined  as 


i9(r)  =  max  0.  min  (2rT  0.25+  0.75r,  0.75  4 ■  0.25r,  2)] . 


QUICK 


(386) 


in  which  the  term  singled  out  by  the  underbrace  is  the  contribution  from  the  QUICK 
scheme  (viz.,  if  i?(r)  were  exactly  equal  to  the  term  highlighted  by  the  underbrace  in  Equa¬ 
tion  (38b),  then  the  advective  flux  approximation  reduces  exactly  to  Leonard’s  QUICK 
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scheme).  Furthermore,  in  Equation  (38a),  the  following  definitions  have  been  used: 


(u)*  =  (u±  |u|)/2.  (u)*  =  (t'±  |v|)/2,  (u>)* 

=  (to±  H)/2; 

(39a) 

and,  r*  ( k  €  {to,  e,  n,  s,  t,  6})  are 

+  _  <t>p  ~  4>w  _  _  4>b  -  4>ee  +  _  <f>w  -  4>ww 

$e  —  4>p  e  «£p  -  4>e  w  4>p  -  4>w 

4>P  - 

r”  ~  4>w-4>p  ’ 

(396) 

+  _  <t>p  —  4>s  _  __  4>n  —  <Pnn  +  _  <f>s  —  <f>ss 

4>n  —  4>p  '  "  <t>p  -  4>n  s  4>p  -  4>s  ' 

"t 
*  1 

II 

^  It 

1  1 

(39c) 

—  0B  -  _<t>T  -  4>TT  +  _  4>B  -  <i>BB 

4>r  ~4>p'  1  4>p  -  4>t  :  b  <t>p~<t>B  ' 

_  4>p  —  4>t 

4>b  —  4>p' 

(39d) 

where  EE  refers  to  the  node  east  of  node  E,  WW  refers  to  the  node  west  of  node  W .  and 
so  forth.  It  is  noted  that  the  higher-order  UMIST  scheme  for  approximation  of  advective 
volume-face  fluxes  uses  a  twelve-point  molecule,  with  the  neighborhood  of  a  cell  center  P 
consisting  of  W,  E,  5,  N,  B,  T,  WW,  EE,  SS,  NN,  BB,  and  TT. 

Pressure-correction  algorithm 

For  an  incompressible  flow,  there  is  no  direct  method  for  specifying  an  equation  for  pres¬ 
sure.  In  this  case,  pressure  needs  to  be  interpreted  as  a  kinematic  variable  {Lagrange 
multiplier)  which  is  determined  indirectly  so  as  to  satisfy  the  continuity  equation  [cf.  Equa¬ 
tion  (4)].  Indeed,  the  lack  of  a  temporal  derivative  in  the  continuity  equation  implies  that 
this  equation  can  be  interpreted  as  a  constraint  for  the  pressure  (rather  than  a  transport  or 
evolution  equation).  In  urbanSTREAM,  the  continuity  equation  is  enforced  indirectly  by 
solving  a  pressure-correction  equation  which,  as  part  of  the  iterative  sequence,  steers  the 
pressure  towards  a  state  in  which  all  mass  residuals  in  the  cells  are  negligibly  small.  The 
two  iterative  schemes  available  in  urbanSTREAM  to  enforce  mass  conservation  through  a 
pressure- correct  ion  algorithm  are  the  Semi-Implicit  Method  for  Pressure- Linked  Equations 
(SIMPLE)  described  in  detail  by  Patankar  [13]  and  the  SIMPLEC  (SIMPLE)- Consistent) 
algorithm  of  Van  Doormal  and  Raithby  [14]. 

The  SIMPLE(C)  algorithm  provides  a  method  for  connecting  the  discretized  forms  of  the 
momentum  and  continuity  equations  to  give  an  equation  linking  the  pressure  correction  at  a 
node  to  its  neighboring  nodes.  In  this  iterative  solution  sequence,  tq  ( i  =  1,  2,  3)  are  initially 
obtained  with  an  estimated  pressure  field.  This  is  continuously  updated  by  reference  to  the 
local  mass  residuals,  which  are  used  to  steer  the  pressure  field  towards  the  correct  level. 
However,  the  fully  collocated  variable  storage  used  in  urbanSTREAM,  in  combination  with 
a  central  differencing  for  pressure,  is  known  to  provoke  checkerboard  oscillations,  reflect  ing 
the  pressure- velocity  decoupling.  More  specifically,  if  a  linear  interpolation  scheme  is  used 
to  determine  the  pressure  on  cell  faces,  the  pressure  gradient  across  a  cell  depends  on 
the  pressure  at  the  nodes  in  the  neighboring  cells  and  is  independent  of  the  pressure  in 
the  current  cell,  leading  to  the  above  mentioned  checkerboard  pattern  in  the  pressure. 
To  avoid  this,  the  widely  used  interpolation  practice  of  Rhie  and  Chow  [15]  is  adopted 
to  interpolate  cell-face  velocities  from  the  nodal  values.  This  method  for  interpolation 
essentially  introduces  a  fourth-order  smoothing  term  based  on  the  pressure,  and  prevents 
the  occurrence  of  spurious  pressure  modes.  The  interpolation  scheme  of  Rhie  and  Chow 
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for  removing  checkerboard  pressure  oscillations  within  a  collocated  scheme  for  solving  the 
RANS  equations  in  a  generalized  curvilinear  coordinate  system  has  been  presented  in  detail 
by  Lien  and  Leschziner  [16J. 


The  key  concept  underlying  Rhie  and  Chow’s  interpolation  for  avoidance  of  checkerboard 
oscillations  in  the  pressure  within  a  collocated  storage  arrangement  can  be  seen  as  follows. 
To  this  purpose,  let  us  consider  the  discretized  u- momentum  transport  equation,  which  by 
virtue  of  Equation  (35a)  [with  the  identification  4>  =  u],  is  given  by 

Snb  ®nb^nt>  4"  Qc  /■_  _  \  Ap  /jn\ 

Up  =  - (pe  ~pw)  - ,  (40) 

Up  Cl  P 

Up  D% 


where  the  discretized  pressure  gradient  term  contribution  to  Qc  has  been  split  out,  and  Ap 
is  the  cross-sectional  area  of  the  cell  at  P.  Equation  (40}  can  be  re-arranged  as  follows: 

up  =  Hp  -  Dp(pe~pw).  (41) 

A  corresponding  expression  can  be  written  for  the  neighboring  node  to  the  east,  so 

=  He  -  DE:{pee  -  Pe),  (42) 

where  ee  denotes  the  eastern  face  of  the  celt  that  has  node  E  at  its  center.  Similarly,  the 
velocity  at  the  eastern  face  of  the  cell  that  has  node  P  at  its  center  is  given  by 

ue  =  He-  D“  (p£  -  pp ) .  (43) 


The  critical  point  to  note  is  that  the  face  velocity  ue  is  evaluated  by  a  linear  interpolation  of 
the  nodal  velocities  up  and  up  from  which  then  the  pressure-gradient  terms  are  subtracted 
and  to  which  a  compensating  pressure-gradient  fragment  is  added,  the  last  formed  only 
with  the  two  nodal  pressures  that  straddle  the  face  velocity  ue.  The  result  gives 

™e  =  l(Hp  +  HB)-±(D%+D%)tpB-pP) 

=  \  up  +  pw) +uE  + D%(pee -pe)}~  ^(D%  + D%)(pE  -  pp),  (44) 

which  can  be  recast  into  the  following  more  perspicuous  form: 


Ue  =  ^  (Up  +  U£) 

C - - - ' 

Linear  interpolation 


1 

+  2 


D%{pe  ~  pw)  +  D%{pee  ~  Pe)  ~  {D%  +  D%)  {pE  -  pP)]  . 


(45) 


pressure  smoothing 


Equation  (45)  implies  that  the  Rhie  and  Chow  interpolation  practice  consists  of  two  parts:  a 
centered  approximation  for  ue  and  a  stabilizing  pressure  smoothing  term  which  is  a  function 
of  pressure  at  the  neighboring  nodes.  Note  that  in  Equation  (45)  depends  on  the  values 
of  pressure  at  four  nodes,  two  on  either  side  of  the  cell  face. 
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The  pressure-correction  equation  now  arises  by  decomposing  the  correct  pressure  field  p 
into  an  approximate  (guessed)  pressure  field  p*  and  a  corrective  pressure  perturbation  p', 
so  that 

p  —  p*  +  pL  (46) 

Similarly,  the  correct  face  velocity  ue  is  decomposed  into  an  estimated  (or,  guessed)  value 
u*  and  a  correction  u^: 

ue  =  u'e  +  fig,  (47) 

where  u'e  fts  d'*  (p'p  —  p'^,),  with  denoting  the  value  of  du  at  cell  face  e  obtained  as  a 
centered  average  of  the  values  of  du  at  the  two  neighboring  nodes  on  either  side  of  the  face 
(viz.,  at  nodes  P  and  E).  The  form  for  du  depends  on  the  pressure- based  algorithm  used. 
In  the  SIMPLE  algorithm,  =  Aefae  =  .  However,  for  the  SIMPLEC  algorithm  which 

provides  a  better  approximation  for  the  corrective  velocity  perturbation,  tf“  =  Aej (ae  -- 
^anH).  Substitution  of  Equation  (47)  along  with  analogous  expressions  for  the  other  cell 
face  velocities  into  the  continuity  equation,  and  integrating  over  the  volume  of  a  cell,  yields 
the  following  discretization  equation  for  p': 

apPp  =  “nbPnb  +  ***,  (48) 

nb 

where 


4  =  4  Am 

aW  — 

(49a) 

4  =  4a, 

(496) 

4  =  d?At, 

4  =  4A, 

(49e) 

and 

ftp  ~  CLyy  -jh  (X^f  “I"  “h  (49d) 

with  the  mass  imbalance  b ?  given  by 

bP  =  -  u'eAe)  +  -  u*j4n)  +  A)-  (49e) 

We  note  that  if  b $  vanishes  identically,  then  the  starred  velocity  field  satisfies  the  continuity 
equation  and  no  pressure  correction  is  needed.  The  term  b ?  can  be  interpreted  as  a  “mass 
source’1  which  the  pressure  corrections  must  annihilate  (through  the  associated  velocity 
corrections)  in  order  for  the  continuity  equation  to  be  exactly  satisfied. 

The  pressure-correction  equation  for  pf  is  susceptible  to  divergence  unless  some  under¬ 
relaxation  is  applied  during  the  iterative  process.  In  view  of  this,  only  a  fraction  of  pf  is 
added  to  p*  in  Equation  (46)  to  give  the  corrected  pressure  p,  so  that 

P  =  P*  +  OfpP'i  (50) 

where  ap  €  (0,  1].  For  the  SIMPLE  algorithm,  ap  is  typically  chosen  to  have  a  value  of  about 
0.7,  However,  it  is  important  to  note  that  under-relaxation  of  the  pressure  is  not  required 
for  the  SIMPLEC  algorithm,  so  in  this  case  ctp  —  1.0.  Similarly,  we  also  under- relax  the 
corrections  for  velocity  (cf.  Equation  (47)  and  analogous  expressions  for  the  other  velocity 
components).  The  under-relaxation  of  the  velocity  does  not  need  to  be  computed  directly, 
but  rather  it  is  more  convenient  (and  elegant)  to  implement  this  under- relaxation  while 
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solving  the  discretized  equations  for  ut  (i  —  1,  2,3).  With  reference  to  Equation  (35a)  with 
<t>  =  Hi ,  the  under-relaxation  of  the  discretized  ur -momentum  equation  takes  the  following 


form: 


+ Cc  + 


nh 


(l-o*)-*- 


U, 


(51) 


where  €  (0,  1]  is  the  under- relaxation  factor  for  S4;  and,  u*  is  the  value  of  Q*  from 
the  previous  iteration.  In  Equation  (51),  the  contribution  of  the  pressure  gradient  term  has 
been  absorbed  in  Q £  and  the  effect  of  the  under- relaxation  of  the  momentum  equation  is  to 
introduce  an  additional  contribution  to  the  source  term.  Its  effect,  therefore,  is  to  increase 
the  diagonal  dominance  of  the  coefficient  matrix  through  the  inclusion  of  this  additional 
source  term.  Finally,  we  note  that  the  pressure-correction  equation  is  also  affected  by 
the  under- relaxation  of  the  fi*- momentum  equations  in  the  sense  that  the  dUt -terms  in  the 
pressure-correction  equation  [cf.  Equation  (49)]  need  to  be  replaced  as  follows: 


d*'  r-  1,2,3. 


(52) 


The  solution  procedure  for  the  velocity  field  in  urbanSTREAM  is  as  follows.  For  each 
time  step,  the  following  procedure  is  undertaken:  (a)  Guess  the  pressure  field  p *  and  solve 
the  mean  momentum  equations  to  obtain  the  provisional  Cartesian  velocity  components  u * 
(i  —  1,2,3);  (b)  With  the  provisional  velocity  u*  ( i  =  i,2,3),  solve  the  pressure-correction 
equation  to  obtain  p';  (c)  Update  the  provisonal  velocities  and  pressure  to  ut  (i  —  1,2,3) 
and  p  using  the  velocity-correction  and  pressure-correction  formulae;  (d)  Solve  all  the  other 
discretized  transport  equations  (e+gM  A:  and  e):  (e)  Treat  the  corrected  pressure  p  as  the  new 
guessed  pressure  p*;  and,  (f)  Return  to  Step  (a)  until  a  converged  solution  is  obtained  at 
which  point  proceed  to  the  next  time  step. 

Within  this  scheme,  the  transport  equations  for  Ui  {i  =  1,2,3),  kh  and  e  and  the  pressure- 
correction  equation  are  solved  sequentially  and  iterated  to  convergence,  defined  by  reference 
to  Ll -residual  norms  for  the  mass  and  momentum  components.  Here,  the  LI -residual  norm 
is  defined  as  the  sum  of  absolute  residuals  over  all  grid  points  of  the  flow  domain.  The 
residual  norm  provides  a  quantitative  measure  of  how  perfectly  the  discretization  equations 
are  satisfied  by  the  current  values  of  the  dependent  variables.  The  LI -residual  norms  for  the 
mass  and  momentum  components  were  normalized  by  the  mass  and  momentum  fluxes  at  the 
inflow  plane.  A  convergent  solution  was  assumed  after  each  normalized  Ll-resjdiial  norm 
decreased  below  0.001  (although  this  level  for  declaration  of  convergence  can  be  chosen  by 
the  user).  In  urbanSTREAM,  the  discretized  equations  are  solved  using  an  iterative  method 
called  the  strongly  implicit  procedure  (SIP)  that  was  proposed  by  Stone  [17],  The  SIP 
method  uses  an  incomplete  lower-upper  (LU)  decomposition  and  was  specifically  designed 
for  solving  algebraic  equations  that  arise  from  the  discretizations  of  partial  differential 
equations  {e,gM  those  that  arise  in  various  CFD  applications).  For  these  types  of  problems, 
Stone’s  SIP  method  generally  converges  in  a  small  number  of  iterations. 
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Boundary  condition  implementation 


The  boundary  conditions  used  in  urbanSTREAM  (flow  solver)  and  urban(A)£U  (urban 
dispersion)  pertain  to  physical  situations  described  by  partial  differential  equations  of  ei¬ 
ther  a  parabolic  or  elliptic  character;  For  the  mean  concentration  C  and  C*  ,  at  ail  solid 
boundaries  (e.g.,  ground  surface,  walls  and  roofs  of  buildings,  etc.)  in  the  computational 
domain,  zero  Neumann  conditions  were  used,  implying  no  flux  of  C  and  C *  across  the  solid 
boundary  (viz.,  QC fdn\B  —  0  and  dC*  fdn\B  =  0  where  n  refers  to  the  direction  normal 
to  the  solid  boundary  B ).  At  the  computational  (flow)  domain  boundaries,  zero  Neuman 
conditions  were  applied  For  C,  but  convective  boundary  conditions  were  used  for  C*  (viz., 
A'efT dC*/dn  +  ii  ■  n C*  —  0  where  =  (A"  +  Kt)  and  n  is  the  unit  outward  normal  vector 
to  the  flow  domain  boundary).2 

The  boundary  conditions  imposed  on  the  flow  field  are  more  complex.  At  the  inflow  bound¬ 
aries  of  the  computational  domain,  all  transported  flow  properties  need  to  be  prescribed 
(i.e.T  Dirichlet  boundary  conditions  apply  here).  With  the  velocity  prescribed,  no  boundary 
conditions  are  needed  for  the  pressure.  In  urbanSTREAM,  the  mean  velocity  components  tit 
(i  —  1,2,3)  and  the  turbulence  kinetic  energy  k  at  the  inflow  boundary  planes  are  provided 
by  urbanGRID  which  imports  four-dimensional  (3-space  dimensions  plus  1-time  dimension) 
meteorological  fields  from  urban  GEM  LAM  and  interpolates  these  field  quantities  on  the 
grid  nodes  of  the  inflow  planes.  Alternatively,  measured  profiles  of  mean  velocity  and  turbu¬ 
lence  kinetic  energy  (e.g.,  obtained  from  sodars,  radiosondes,  radars,  etc.)  can  also  be  used 
by  urbanSTREAM  to  define  the  boundary  conditions  for  mean  flow  and  turbulence  on  the 
inflow  boundary  planes,  if  these  are  available.  Unfortunately,  the  viscous  dissipation  rate  e 
is  not  available  from  urban  GEM  LAM  and  is  virtually  never  available  from  experimental 
measurements.  In  urbanSTREAM,  e  is  prescribed  on  the  inflow  boundary  planes  using  the 
assumption  of  local  equilibrium  of  the  turbulent  flow  (so  that  the  rates  of  production  and 
destruction  of  turbulence  are  in  near  balance)  to  give 


^3/4  £3/ 2 
min(fcvz,  h/3)  1 

where  z  is  vertical  distance  from  the  ground  and  h  is  the  atmospheric  boundary- layer 
height3. 

Far  downstream,  at  the  outflow  boundary  planes,  the  flow  was  assumed  to  reach  a  fully- 
developed  state  where  no  changes  occur  in  the  flow  direction  (approximately  or  better). 
Hence,  at  the  outflow  boundary,  the  horizontal  gradients  of  all  flow  variables  are  assumed 
to  be  zero;  viz.,  du/dxi  =  dv/dxi  —  dw/dx{  —  dk/dx{  =  defdxi  —  0  (i  =  1,2).  At  the 
upper  boundary,  we  used  free-slip  conditions  for  all  flow  variables.  Furthermore,  at  the 
outflow  boundaries,  a  mass  flux  correction  is  undertaken  to  ensure  that  the  total  mass  flux 
that  exits  through  the  outflow  boundary  planes  is  equal  to  the  total  mass  flux  that  enters 
through  the  inflow  boundary  planes  (and,  in  fact,  the  imposition  of  this  condition  enhances 
convergence  of  the  iterative  flow  solver).  During  iterations  of  the  pressure- based  algorithm 

The  convective  boundary  conditions  need  to  be  used  for  C*  to  ensure  that  C *  is  the  adjoint 
of  C,  verifying  the  duality  relationship  of  Equation  (31). 

The  upper  limit  of  hj 3  in  the  turbulence  length  scale  used  in  Equation  (53)  is  consistent 
with  that  used  in  the  limited- length-scale  fe-£  turbulence  model. 
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(either  SIMPLE  or  SIMPLEC),  there  is  no  guarantee  that,  the  exit  velocities  through  the 
outflow  boundary  planes  will  conserve  mass  over  the  computational  domain  as  a  whole.  To 
ensure  that  the  continuity  is  satisfied  globally,  the  total  mass  flux  M*’  entering  into  the 
domain  and  total  mass  flux  M*;it  going  out  of  the  domain  in  the  x^-coordinate  direction 
(i  =  1,2)  are  first  calculated  as  follows: 


Mi. 


J€l*‘  jeO*‘ 


*  =  1,2, 


(54) 


where  TXi  and  0Xt  denote  the  inflow  and  outflow  boundary  planes  normal  to  the  x,- 
coordinate  direction,  respectively;  and,  (u,)^  corresponds  to  the  inlet  (exit)  velocity  through 
inlet  (exit)  cell  j  whose  inlet  (exit)  face  area  is  A Aj.  There  is  no  unique  procedure  for  en¬ 
suring  that  the  mass  flux  out  is  equal  to  mass  flux  in  M*n'  for  i  —  1,2.  One  method 
to  ensure  overall  continuity  is  to  add  a  constant  velocity  correction  ti'  given  by 


u 


C 

i 


(*C  -  AC,) 


(55) 


to  the  outlet  velocity  component  j  E  Ox'  ( i  =  1,2),  Alternatively,  global  continuity 

in  the  computational  domain  can  also  be  ensured  by  multiplying  the  outflow  plane  velocities 
{&i)j  (j  €  Qx')  by  the  ratio  Currently,  urbanSTREAM  uses  the  first  method 

for  ensuring  that  mass  is  conserved  over  the  entire  computational  domain. 


At  all  the  solid  boundaries  (ground,  obstacle  walls,  obstacle  roofs),  standard  wall  functions 
are  applied  for  the  mean  velocities  and  turbulence  quantities.  Wall  functions  are  used 
to  reduce  the  computational  cost  associated  with  the  alternative  of  using  a  low- Reynolds 
number  turbulence  model  to  numerically  integrate  the  solution  through  the  entire  near- 
wall  region,  including  the  viscous  sublayer  up  to  the  wall  where  no-slip  and  impermeability 
conditions  can  be  applied.  Typically,  this  will  require  a  very  fine  grid  near  the  wall  because 
spatial  variations  in  near- wall  turbulence  structure  are  large  here  due  to  the  combined 
influence  of  viscosity  and  wall-induced  anisotropy.  Log- law  based  4 wall  laws'  are  adopted  in 
conjunction  with  high- Reynolds  number  turbulence  models  to  bridge  the  viscous  sublayer. 
Basically,  the  flow  quantities  at  the  first  grid  point  above  the  solid  surface,  which  is  located 
outside  the  viscous  sublayer,  are  related  to  the  wall  friction  velocity  based  on  the  assumption 
of  a  semi-logarithmic  velocity  distribution. 


The  implementation  of  a  wall  function  boundary  condition  for  the  mean  velocity  tangential 
to  the  wrall  first  requires  the  evaluation  of  the  normalized  distance  perpendicular  to  the  wall: 


z 


+ 


(56) 


where  zp  is  the  distance  from  the  near- wall  node  P  to  the  solid  surface,  rw  is  the  wrall  shear 
stress,  and  p  is  fluid  density*  If  z+  <  11.6,  the  flow  at  near-wall  node  P  is  assumed  to  be  in 
the  laminar  sublayer  where  the  flow  is  linear:  otherwise,  the  flow  is  turbulent,  lying  in  the 
log-law  region  where  a  log  variation  in  the  tangential  velocity  is  assumed  to  prevail  and  a 
wall  function  is  used.  For  the  velocity  component  *13  normal  to  the  wall*  the  impermeability 
condition  that  U3  —  0  at  the  wall  is  used* 
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More  specifically,  if  z+  <  11.6.  the  wall  shear  stress  in  the  x.-coordinate  direction  (i  =  1, 2) 
is  obtained  from  r*1  =  pu(ui)p/zp ,  where  (u,)p  is  the  tangential  velocity  at  the  near¬ 
wall  node  P  (assuming  a  linear  variation  of  velocity  with  distance  from  the  wall).  The 
shear  force  Fa  on  the  bottom  face  of  the  near-wall  cell  (with  center  node  P)  is  given  by 
Fs  =  —pv{ux)pAb/zp  (,4b  is  the  wall  (bottom)  area  of  the  near-wall  cell).  The  link  to  the 
wall  in  the  discretized  Uj-momentum  equation  (r  =  1,2)  is  suppressed  by  setting  ap  =  0 
and  the  wall  function  contribution  S'p  =  —puAb/zp  is  added  to  Sp  (cf.  Equation  (37)]. 
If,  however,  z+  >  11.6,  the  near-wall  node  P  is  considered  to  lie  in  the  log-law  region  of 
the  turbulent  flow  and  a  wall  function  is  utilized.  To  this  purpose,  one  calculates  the  wall 
shear  stress  tw  and  the  associated  wall  shear  force  Fs  from  the  log-law  velocity  profile  with 
properties  evaluated  at  the  near-wall  node  P.  Again,  the  coefficient  ap  in  the  discretized 
u.-momenturn  equation  (i  =  1,2)  is  nullified  and  the  source  Qc  [cf.  Equation  (35b)]  is 
modified  to  explicitly  include  the  contribution  from  the  shear  force  acting  on  the  bottom 
near-wall  cell  face  as  follows: 


Qc 


Qc- 


pk^Cl/% 

\n(Eklp2  zp  / 1/) 


(^i)  pAb ' 


i  =  1.2, 


(57) 


where  E  =  9.793  is  a  constant. 

To  implement  the  wall  function  of  Equation  (57),  a  correct  value  for  kp  (TKE  at  the 
near-wall  node)  is  required.  A  local  equilibrium  of  turbulence  (where  the  volume-averaged 
production  and  dissipation  of  turbulence  are  in  balance)  is  assumed  to  prevail  near  the  wall. 
The  average  production  of  A;  in  a  near-wall  cell  is  given  by 


Pk  =  \Tw\-\up\/zP,  (58) 

where  Ir^l  =  ((r*1)2  +  (r*1)2)'72  and  |u£P|  =  ((ui)P  +  (u2)p)'/2  are  the  magnitudes  of 
the  wall  shear  stress  and  the  tangential  velocity  (at  near- wall  node  P).  respectively.  With 
respect  to  implementation,  the  net  fc-source  per  unit  volume  in  Equation  (58)  is  included  in 
the  near-wall  cells  in  the  discretized  fc-transport  equation  by  suppressing  the  link  to  the  wall 
boundary  (ap  —  0)  and  augmenting  the  source  term  Qc  with  the  additional  contribution 
gpr°d  _  |T^|  .  (All  is  the  volume  of  the  near-wall  cell).  Finally,  the  viscous 

dissipation  e  at  a  near-wall  cell  with  center  node  P  has  the  nodal  value  set  at  (with  an 
assumed  local  equilibrium  condition  persisting  in  the  near-wall  cell) 

eP  =  C3J4k3J2/(kvzP).  (59) 


Application  to  real  urban  environment 


Here,  a  computational  study  is  performed  in  which  the  capabilities  of  urbanSTREAM  for 
urban  flow  prediction  and  urban(A)EU  for  urban  dispersion  predictions  are  examined  by  ref¬ 
erence  to  a  full-scale  experimental  study  of  flow  and  dispersion  in  a  real  cityscape  (Oklahoma 
City,  Oklahoma). 
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Joint  Urban  2003 


Joint  Urban  2003  (JU2003)  experiment  was  a  major  urban  study  that  was  conducted  in 
Oklahoma  City,  Oklahoma  during  the  period  from  June  28  to  July  31,  2003  [18].  The 
principal  objective  of  JU2003  was  to  obtain  high-quality  meteorological  and  tracer  data 
sets  documenting  urban  flow  and  dispersion  in  a  real  city  on  a  range  of  scales;  namely,  from 
flow  and  dispersion  in  and  around  a  single  city  block  (street  canyon  at  building  scale),  to  that 
in  and  around  several  blocks  in  the  central  business  district  (CBD)  of  downtown  Oklahoma 
City  (neigborhood  scale),  and  finally  to  that  in  the  suburban  area  several  kilometers  from 
downtown  Oklahoma  City  (urban  scale).  Additionally,  JU2003  included  measurements 
of  indoor  flow  and  dispersion  in  four  buildings  in  the  CBD  that  were  coordinated  and 
conducted  in  conjunction  with  the  outdoor  urban  field  experiments  in  order  to  study  physical 
mechansims  involved  in  indoor-outdoor  exchange  rates. 

A  large  number  of  meteorological  measurements  in  Oklahoma  City  were  obtained  during 
JU2003.  These  included  detailed  measurements  of  wind  and  turbulence  characteristics  in 
the  urban  atmospheric  boundary  layer  obtained  with  remote  sensing  instruments  such  as 
Doppler  lidars,  Doppler  sodars  (acoustic  wind  profiler)  and  radar  profilers;  and,  measure¬ 
ments  of  wind,  temperature,  and  turbulence  in  the  urban  canopy  layer  (including,  within  a 
street  canyon)  using  fast-response  in-situ  meteorological  sensors  such  as  sonic  anemometers 
and  infrared  thermometers.  Furthermore,  tracer  bag  samplers  were  used  to  measure  mean 
concentration  data  obtained  from  the  release  of  a  sulfur  hexafluoride  (SFg)  tracer  in  down¬ 
town  Oklahoma  City,  Oklahoma  at  three  different  locations.  To  this  purpose,  ten  intensive 
observation  periods  (IOPs)  were  undertaken  during  JU2003,  during  each  of  which  there 
was  typically  three  30-minute  (continuous)  tracer  gas  releases  as  well  as  four  puff  (instanta¬ 
neous)  releases  where  balloons  filled  with  tracer  gas  were  “popped".  For  all  the  IOPs,  eacli 
of  the  three  30-minute  continuous  releases  was  separated  by  90  minutes  and  this  set  of  three 
continuous  releases  was  either  preceded  or  followed  by  a  set  of  four  instantaneous  releases 
with  a  20-minute  interval  between  each  release.  The  release  rates  for  the  continuous  source 
wets  constant  (generally  to  within  about  5%)  during  each  experiment,  and  the  release  rates 
from  the  various  experiments  ranged  over  the  interval  from  2--5  g  s-1.  The  tracer  masses 
released  from  the  instantaneous  source  varied  over  the  range  from  300-1000  g.  The  tracer 
gas  from  all  releases  was  sampled  in  and  around  downtown  Oklahoma  City  on  the  regular 
CBD  sampler  grid  and  as  far  downwind  as  four  kilometers  from  the  release  along  various 
sampling  arcs.  The  first  six  IOPs  were  conducted  during  the  day  from  08:00  CST  to  16:00 
CST.  whereas  the  last  four  IOPs  were  conducted  during  the  night  from  22:00  CST  to  06:00 
CST. 

Grid  generation  for  Oklahoma  City 

Some  capabilities  of  the  modeling  system  are  presented  with  the  aid  of  numerical  simulations 
performed  for  a  built-up  (urban)  region  in  Oklahoma  City,  Oklahoma.  The  modeling  domain 
with  the  extent  of  1,  934,25  m  x  3,  610,6  m  x  800.0  m  in  the  x-  (or,  W-E),  y-  (or.  S-N)  and 
z-  (or,  vertical)  directions,  respectively,  covers  the  CBD  of  Oklahoma  City  and  surrounding 
environs.  At  ground  level  where  z  —  0,  the  southwest  corner  of  the  modeling  domain 
(see  Figure  6)  is  at  the  following  coordinates  in  the  Universal  Transverse  Mercator  (UTM) 
coordinate  system4:  zone  =  14.  .To  —  633.683  UTM  easting  and  yo  =  3,923.940  UTM 

The  UTM  easting  coordinate  reported  here  and  elsewhere  in  this  report  is  referenced  relative 
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northing  (or,  equivalently,  in  the  geodetic  coordinate  system  this  location  is  35.449959°  N 
and  -97.52694°  E).  The  internal  coordinate  system  used  in  urbanSTREAM  is  shown  in 
Figure  6,  where  the  southwest  corner  of  the  modeling  region  is  chosen  as  the  origin  (0, 0)  in 
the  i -y  (horizontal)  plane.  All  distances  shown  tiere  have  been  normalized  by  a  reference 
length  scale  which  is  chosen  in  this  case  to  be  Aref  =  644.75  m.  Hence,  in  this  internal 
coordinate  system,  the  northeast  corner  of  the  modeling  region  is  referenced  as  (3,5,6).  A 
proper  subset  within  this  modeling  region  is  chosen  as  the  region  in  which  buildings  will  be 
explicitly  resolved  in  the  flow  simulation;  for  this  example,  this  rectangular  building-aware 
region  (644.75  m  x  709.23  m)  has  its  southwest  corner  at  (1. 2.5)  and  its  northeast  corner  at 
(2,3.6).  In  the  portion  of  the  modeling  region  lying  outside  the  building-aware  region,  all 
buildings  are  treated  as  virtual  and  their  effects  on  the  flow  are  modeled  using  a  distributed 
drag  force  representation  in  the  mean  momentum  equations. 

ESRI  Shapefiles  of  the  shapes,  locations  and  heights  of  buildings  in  Oklahoma  City  are 
available  from  the  Joint  Urban  2003  archival  database  (https://ju2003-dpd.dpg.army.inil). 
These  Shapefiles  were  used  in  urbanGRID  to  generate  automatically  a  grid  mesh  over  the 
modeling  region  as  shown  in  Figure  7.  The  simulation  was  carried  out  in  a  three-dimensional 
Cartesian  framework,  and  curved  surfaces  on  buildings  or  planar  building  surfaces  that 
are  not  aligned  with  the  grid  tines  are  approximated  by  stepwise  surfaces.  A  mesh  of 
99  x  139  x  69  grid  lines  in  the  x-,  y-,  and  z-directions,  respectively,  was  used  to  accomodate 
all  the  necessary  geometrical  details.  The  interior  building- aware  region  was  covered  with  a 
fine  calculation  grid  of  55  x  100  x  69  grid  lines  to  better  approximate  the  building  features 
in  this  region.  The  grid  arrangement  adopted  here  is  shown  in  the  x-y  plane  in  Figure  7. 
Hence,  the  fine  grid  used  for  the  building- aware  region  contains  379,500  nodes,  whereas 
the  entire  computational  domain  was  covered  with  a  mesh  of  945,509  nodes.  The  grid 
lines  were  preferentially  concentrated  near  the  solid  surfaces  (ground,  building  rooftops  and 
walls)  where  the  gradients  in  the  flow  properties  are  expected  to  be  greatest,  and  the  spacing 
between  the  grid  lines  was  gently  stretched  with  increasing  distance  from  the  solid  surfaces. 
Figure  7  shows  the  finite-volume  approximation  of  the  explicitly  resolved  buildings  in  the 
building-aware  region,  and  these  approximated  buildings  should  be  compared  to  the  actual 
(true)  buildings  in  this  area  given  in  Figure  8  in  the  form  of  a  satellite  photograph  and  a 
three-dimensional  computer  rendering  of  the  buildings  from  the  ESRI  Shapefile.  As  can  be 
seen,  there  is  very  little  difference  between  finite  volume  approximated  buildings  and  the 
actual  (true)  buildings  and  therefore  the  stepwise  approximation  of  the  building  surfaces  is 
not  expected  to  undermine  the  accuracy  of  the  subsequent  flow  simulation. 

Flow  field 

The  flow  field  in  the  computational  domain  was  computed  using  urbanSTREAM  in  a  stand¬ 
alone  mode  (viz.,  this  flow  model  was  not  coupled  to  the  urban  GEM  LAM  model).  In 
consequence,  at  the  inflow  (inlet)  boundary  of  the  computational  domain,  the  measured 
profiles  of  the  undisturbed  mean  velocity  and  turbulence  kinetic  energy  are  used.  The 
flow  simulation  conducted  here  was  for  IOP-9  for  the  time  period  from  06:00-06:30  UTC 
(01:00-01:30  CDT)  on  28  July  2003  when  the  prevailing  winds  were  from  the  south  at 
about  6.8  m  s'1  at  50-rn  above  ground  level  at  the  southern  edge  of  the  computational 
domain  shown  in  Figure  6.  As  a  result,  the  inflow  mean  velocity  and  turbulence  kinetic 
energy  profiles  (required  for  the  inflow  specification)  were  obtained  from  the  Scintec  Multi- 

to  the  central  meridian  of  the  zone. 
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frequency  Flat  Array  Antenna  Series  (MFAS)  soclar  operated  by  a  team  of  investigators  from 
Pacific  Northwest  National  Laboratory  (PNNL).  The  Scintec  MFAS  sodar  has  an  operating 
frequency  in  the  2  kHz  band  allowing  for  a  spatial  resolution  of  10  m  with  a  measurement 
range  of  up  to  1000  m.  This  sodar  was  located  2  krn  south  of  the  central  business  district 
of  Oklahoma  City  in  a  parking  lot  at  the  Traffic  Maintenance  Yard  (latitude  35.45°  N: 
longitude  —97.53°  E)  which  is  on  the  south  edge  of  Wheeler  Park.  The  area  south  of 
this  sodar  was  covered  with  low-rise  (less  than  three  stories)  buildings,  houses,  and  trees. 
The  PNNL  sodar  was  located  very  near  the  upwind  (southern)  edge  of  the  computational 
domain,  and  measurements  of  mean  wind  velocity  and  turbulence  were  used  to  define  the 
inflow  conditions  for  urbanSTREAM. 

The  flow  simulation  in  the  computational  domain  was  carried  out  using  urbanSTREAM  in 
the  unsteady  RANS  mode.  The  effects  of  the  buildings  outside  the  building- aware  region 
were  represented  by  a  distributed  drag  force  approximation  in  the  mean  momentum  equa¬ 
tion  with  a  normalized  drag  coefficient  of  Cd  =  CdAAref  =  100  (where  Aref  is  the  reference 
length  scale).  For  the  simulations,  a  constant  time  step  At  =  30  s  was  used.  The  predicted 
wind  statistics  were  averaged  over  500  time  steps  after  the  flow  had  achieved  a  (pseudo) 
steady  state.  We  compared  the  predicted  vertical  profiles  of  the  mean  horizontal  wind  speed 
with  associated  measured  values  obtained  at  three  different  locations  in  the  computational 
domain.  These  comparisons  are  shown  in  Figure  9.  The  measurements  of  wind  speed 
were  obtained  from  three  different  sodars:  two  minisodars  deployed  by  Argonne  National 
Laboratory  (ANL)  and  one  sodar  operated  by  National  Atmospheric  and  Oceanic  Admin¬ 
istration  (NOAA)  Atmospheric  Research  Laboratory  Field  Research  Division  (ARLFRD). 
One  ANL  minisodar  was  located  near  the  southern  edge  of  the  central  business  district  in 
the  Oklahoma  City  Botanical  Garden  (latitude  35.46°  N;  longitude  —97.52°  E).  The  CBD 
of  Oklahoma  City,  which  includes  many  high-rise  buildings,  was  immediately  north  of  this 
sodar.  The  second  ANL  sodar  was  located  approximately  1  km  north  of  the  CBD  of  Okla¬ 
homa  City  near  the  Minor  building  (latitude  35.48°  N;  longitude  —97.52°  E).  The  ARLFRD 
sodar  (Radian  600PA  phased-array  Doppler  sodar)  was  located  about  1.5  km  northeast  of 
the  CBD  of  Oklahoma  City  at  the  campus  of  the  Oklahoma  School  of  Science  and  Math¬ 
ematics  on  the  southwest  corner  of  Lincoln  and  NE  13th  Street  (latitude  35.48147°  N  ; 
longitude  —97.5051°  E).  Both  ANL  minisodars  were  operated  with  the  first  range  gate  at 
5  tn,  a  range  gate  spacing  of  5  m  with  the  top  gate  at  200-m  above  ground  level  (ACL), 
and  an  averaging  (integration)  time  of  15  min.  The  ARLFRD  sodar  was  operated  with  the 
first  range  gate  at  40  m,  a  range  gate  spacing  of  10  m  with  the  top  gate  at  300-m  AGL, 
and  an  averaging  (integration)  time  of  15  min. 


From  Figure  9,  it  is  seen  that  the  model  with  the  inclusion  of  the  drag  force  representation 
for  virtual  buildings  (dashed  curves)  correctly  reproduces  the  overall  qualitative  trends  in 
the  horizontal  wind  speeds  at  the  three  locations  in  the  computational  domain.  At  the 
locations  of  the  ANL  minisodars.  the  wind  speeds  predicted  by  the  model  were  slightly 
larger  than  the  measured  wind  speeds  over  most  of  the  range  between  about  25  m  and 
200  m  AGL.  The  predicted  wind  speeds  agreed  well  with  those  measured  at  the  location  of 
the  ARLFRD  sodar,  although  the  measurements  here  from  one  15-min  averaging  period  to 
another  showed  greater  variability  than  those  measured  at  the  two  ANL  minisodar  locations. 
Figure  9  also  shows  model  predictions  of  horizontal  wind  speed  without  the  inclusion  of  a 
distributed  drag  force  approximation  for  the  virtual  buildings  (solid  lines).  Note  in  this  case 
that  the  wind  speed  is  greatly  overestimated  by  the  model  owing  to  the  fact  that  without 
a  distributed  mean-momentum  sink  within  the  urban  canopy  to  represent  the  effects  of  the 
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unresolved  buildings  on  the  flow,  there  cannot  be  a  deceleration  of  the  flow  in  this  region 
leading  to  the  overprediction  of  the  wind  speed. 


Urban  dispersion 

The  flow  field  statistics  predicted  by  urbanSTREAM  were  next  used  to  “drive1  an  urban 
dispersion  model  in  both  the  source-oriented  (urbanEU)  and  receptor-oriented  (urban AEU) 
modes  within  the  Eulerian  framework.  Furthermore,  these  wind  field  statistics  were  also 
used  to  “drive*'  a  first-order  Lagrangian  Stochastic  (LS)  model  for  urban  dispersion  (ur- 
banLS,  or  more  precisely  urbanLS-1).  The  simulations  were  conducted  for  the  second 
continuous  30-mi  n  release  of  SFq  in  I  OP- 9,  which  occurred  in  the  period  from  06:00-06:30 
UTC  (01:00-01:30  CDT)  on  28  July  2003.  The  SF$  tracer  gas  was  released  at  a  point  from 
a  specially  designed  gas  dissemination  system  that  included  a  mass  flow  controller  attached 
to  a  Campbell  Scientific  CR-23  data  logger  [19].  The  dissemination  point  was  located  on 
the  south  side  of  Park  Avenue  (latitude  35.4687°  N;  longitude  —97.5156°  E)  with  a  release 
height  of  L9  m,  This  source  location  was  near  the  center  of  the  computational  domain  (see 
Figure  6),  The  constant  gas  release  rate  for  this  experiment  was  2.0  gs"1. 

The  SFs  plume  concentrations  were  measured  using  Programmable  Integrating  Gas  Sam¬ 
plers  (PIGS)  deployed  by  ARLFRD.  The  subsequent  analysis  of  the  samples  was  performed 
using  an  Automated  Tracer  Gas  Analysis  System  ( ATG  AS)  which  used  gas  chromatography 
(GC)  analysis  techniques  along  with  autosampler  capabilities  to  give  time- integrated  con¬ 
centration  measurements.  The  sampling  system  was  designed  to  provide  average  SF6  plume 
concentrations  over  specific  time  intervals  at  given  receptor  locations.  The  PIGS  collected 
12  samples  by  pumping  air  into  12  individual  Tedlar  bags.  Subsequently,  the  bag  samples 
were  analyzed  using  the  ATG  AS,  The  grid  of  sampling  stations  at  which  predictions  of  SF^ 
plume  concentrations  were  compared  with  measurements  is  shown  in  Figure  10  in  relation¬ 
ship  to  the  location  of  the  source.  The  samplers  used  for  comparison  include  19  samplers 
located  on  the  CBD  sampling  grid  and  six  samplers  located  along  the  1-km  sampling  arc. 

Figure  1 1  displays  various  isopleths  (on  a  logarithmic  scale)  of  the  predicted  mean  con¬ 
centration  field  C  at  ground  level  in  the  computational  domain  obtained  using  urban EU 
(source-oriented  approach).  Figure  12  shows  isopleths  (on  a  logarithmic  scale)  of  the  influ¬ 
ence  function  C*  (viz.,  adjunct  of  the  concentration  C)  at  ground  level  corresponding  to  one 
of  the  sampling  locations  (receptors)  along  the  1-km  sampling  arc.  This  result  was  obtained 
using  urban  AEU  (receptor-oriented  approach).  The  influence  function  C*  characterizes  the 
concentration  “seen”  at  the  given  receptor  through  the  duality  relationship  of  Equation  (31) 
for  an  arbitrary  source  distribution  Q.  In  essence,  the  value  of  C*  at  a  given  spatial  position 
provides  information  on  the  contribution  of  a  point  source  with  a  unit  emission  rate  to  the 
concentration  at  the  receptor. 


Figure  13  compares  predictions  of  the  mean  concentration  [in  parts-per-trillion  by  volume 
(pptv)|  obtained  using  urbanEU,  urban  AEU  and  urbanLS-1  with  the  experimental  concen¬ 
tration  data  measured  at  10  different  sampling  locations  along  Kerr  Avenue  and  McGee 
Avenue.  Similarly,  Figure  14  shows  comparisons  of  mean  concentration  obtained  using  the 
three  dispersion  models  with  experimental  concentration  measurements  made  at  eight  sam¬ 
pling  stations  located  along  4th  Street  and  5th  Street,  Finally,  Figure  15  exhibits  predicted 
and  observed  mean  concentration  obtained  at  two  sampling  locations  along  6th  Street  and 
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six  sampling  stations  along  the  1-krn  sampling  arc.  The  experimental  concentration  data 
shown  here  is  for  a  30-min  averaging  time.  Generally  speaking,  the  predictions  for  mean  con¬ 
centration  at  or  near  the  mean  plume  centerline  were  quite  good,  with  predictions  within 
a  factor  of  two  of  the  observed  concentration.  However ,  the  predicted  concentrations  at 
sampling  locations  56,  66,  and  76  were  more  than  a  factor  of  five  lower  than  the  experimen¬ 
tally  measured  values,  suggesting  that  either  the  predicted  plume  was  too  narrow  or  the 
eastern  edge  of  the  predicted  plume  was  too  far  west  at  these  locations.  Furthermore,  it  is 
interesting  to  note  that  the  30-min  average  concentration  at  the  receptor  calculated  using 
the  influence  function  {receptor-oriented  approach)  generally  agrees  well  with  the  concen¬ 
tration  predicted  using  the  source-oriented  approach.  The  discrepancy  in  the  predictions 
using  these  two  approaches  is  due  to  the  non-uniform  grid  utilized  in  the  simulations  (the 
mesh  was  finer  in  the  region  surrounding  the  source  and  coarser  generally  in  the  region 
downwind  of  the  source  where  the  receptors  were  located). 


Conclusions 


This  report  provides  a  technical  description  of  the  models  that  comprise  Component  1  of 
CRT  I  Project  02-0093  RD  whose  principal  objective  is  the  development  of  an  advanced, 
fully  validated,  state-of-the-science  modeling  system  for  the  prediction  of  urban  Row  (i,e., 
turbulent  flow  through  cities)  and  the  concomitant  problem  of  modeling  the  dispersion 
of  CBRN  agents  released  in  a  populated  urban  complex.  Component  l  focuses  on  the 
development  of  an  urban  microscale  building- aware  flow  model  and  Eulerian- based  urban 
dispersion  models.  This  capability  has  been  developed  because  almost  all  current  dispersion 
models  employ  simple  diagnostic  wind  fields  and  Gaussian  diffusion  techniques  to  predict 
plume  (cloud)  dispersion  which  are  not  applicable  in  the  highly  disturbed  Hows  present  in 
an  urban  environment. 

The  principal  module  of  Component  1  is  urbanSTREAM,  which  is  a  general  second-order 
accurate  finite- volume  code  designed  for  the  simulation  of  urban  flow  using  a  two-equation 
turbulence  closure  model  (namely,  the  standard  Ar-e  model  and  limited- length-scale  k~e 
model).  This  is  perhaps  the  simplest  complete  turbulence  model  (in  the  sense  that  no 
advance  knowledge  of  any  property  of  the  turbulence  is  required  for  the  simulation  other 
than  the  initial  and/or  boundary  conditions  for  the  problem)  that  is  available  currently  as  a 
general  purpose  simulator  for  urban  flows.  The  modeling  scheme  used  here  is  simple  enough 
to  be  tractable  numerically  and,  hence,  not  require  excessive  computing  time  which  is  im¬ 
portant  if  the  system  is  to  be  employed  for  emergency  response  applications.  In  addition, 
Component  1  incorporates  a  module  (urbanGRID)  for  the  automatic  generation  of  grids 
in  the  computational  domain  when  provided  with  detailed  geometric  information  on  the 
shapes  and  locations  of  buildings  in  the  urban  environment  in  the  form  of  ESR1  Shapefiles. 
Finally,  Component  1  also  includes  modules  for  the  prediction  of  urban  dispersion  in  the 
Eulerian  framework:  namely,  urbanEU  which  is  an  Eulerian  grid  dispersion  model  based  on 
numerical  solution  of  a  K- theory  advection-diffusion  equation  (source-oriented  approach) 
and  urbanAEU  which  is  a  receptor-oriented  dispersion  model  based  on  numerical  solution 
of  the  adjoint  of  the  AT- theory  advect ion-diffusion  equation. 

The  microscale  urban  flow  and  dispersion  modeling  system  of  Component  1  has  been  ap¬ 
plied  in  stand-alone  mode  to  simulate  flow  and  SFg  tracer  dispersion  in  Oklahoma  City, 
Oklahoma.  The  numerical  simulations  of  IQP-9  in  the  Joint  Urban  2003  experiment  pro- 
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vide  an  initial  demonstration  that  the  developed  modeling  system  can  correctly  reproduce 
many  features  of  the  flow  and  dispersion  in  a  real  cityscape.  While  the  results  from  this 
preliminary  case  study  are  certainly  not  definitive*  they  do  suggest  that  the  urban  disper¬ 
sion  modeling  system  has  the  potential  to  improve  the  predictive  performance  capabilities 
for  emergency  response  in  built-up  areas  where  the  flow  is  highly  disturbed. 

The  numerical  models  included  in  Component  1  require  further  improvements  and  eval¬ 
uation  against  data  from  full-scale  meteorological  and  dispersion  field  experiments  in  the 
urban  environment.  In  particular*  an  evaluation  of  the  predictive  performance  of  the  urban 
microscale  flow  and  dispersion  models  when  fully  coupled  with  urban  GEM  LAM  will  be 
required.  An  improvement  of  the  modeling  system  developed  under  Component  1  can  be 
expected  from:  (a)  the  inclusion  of  more  sophisticated  turbulence  closure  schemes  (e+g., 
second- moment  closure  which  potentially  can  reproduce  effects  of  streamline  curvature,  ro¬ 
tation  and  swirl,  secondary  motion*  and  other  effects  in  highly  disturbed  complex  flows 
better  than  with  an  eddy- viscosity  concept);  (b)  application  of  hybrid  Reynolds- averaged 
Navier-Stokes/ large-eddy  simulation  (RANS/LES)  approach  which  allows  very  large,  coher¬ 
ent,  or  deterministic  structures  in  the  urban  flow  to  be  resolved  rather  than  modeled;  (c) 
inclusion  of  thermal  effects  on  flow  and  contaminant  dispersion  in  the  urban  street  canyons; 
and,  (d)  inclusion  of  an  adaptive  mesh  refinement  capability  in  the  models  that  would  allow 
a  dynamic  tracking  of  the  urban  flow  features  as  the  computation  proceeds,  giving  complete 
control  of  grid  resolution  with  an  expected  concomitant  computational  savings  over  a  static 
grid  approach.  However,  we  emphasize  that  these  future  developments  must  adhere  to  the 
underlying  requirement  for  maintaining  a  well-balanced  compromise  between  sophistication 
of  mathematical  modeling  and  availablity  of  computational  resources. 
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Figure  1.  Relationship  between  various  components  of  CRTI  Project  02-0093RD. 
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Figure  2.  Various  modules  of  Component  I  and  their  relationship  to  other  components  ofCRTI  Project  02-0093RD. 
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Figure  3.  Illustration  of  the  ray-casting  approach  to  determine  if  a  grid  cell  with  centroid  c 
lies  inside  the  building  B . 
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Figure  4.  Determination  of  whether  two  line  segments  intersect  each  other. 
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Figure  5.  Finite  volume  and  storage  arrangement. 
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Figure  6.  Computational  domain  used  for  simulation  of  flow  in  Oklahoma  City. 
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Figure  7.  Computational  grid  generated  by  urbanGRID  for  Oklahoma  City. 
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Figure  9.  Comparison  of  mean  wind  speed  measured  at  three  locations  in  the  computational  domain  with  model 
predictions.  The  inflow  conditions  for  flow  model  were  obtained  from  the  measurements  made  by  PNNL  sodar. 
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Figure  10.  Detectors  (squares)  were  positioned  in  the  CBD  sampling  grid  and  at  the  1-km  sampling  arc  for 
measurement  of  the  time-averaged  tracer  concentrations  released  from  the  indicated  source  (solid  circle). 
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side  of  Park  Avenue  (Oklahoma  City). 


the  1-km  sampling  arc  (Oklahoma  City). 
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Figure  13.  Comparison  of  predicted  mean  concentrations  obtained  from  urbanEU,  urban  A  EU  and  urbanLS-1 
with  experimental  measurements  obtained  with  detectors  along  Ken  Avenue  and  McGee  Avenue. 
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Figure  14.  Comparison  of  predicted  mean  concentrations  obtained  from  urbanEU,  urbanAEU  and  urbanLS-1 
with  experimental  measurements  obtained  with  detectors  along  4th  Street  and  5lh  Street. 
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Figure  15.  Comparison  of  predicted  mean  concentrations  obtained  from  urbanEU,  urbanAEU  and  urbanLS-1 
with  experimental  measurements  obtained  with  detectors  along  6,h  Street  and  the  1-km  sampling  arc. 
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